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PREFACE 

Th.e  following  report  has  been  prepared  in  accordance  with  the 
terms  of  Contract  No.  DA-36-03^-OED-1330  and  constitutes  the  Final  Report 
called  for  under  the  terms  of  that  contract.   It  is  divided  into  three 
parts  covering  the  salient  engineering  work  as  well  as  the  mathematical 
and  meteorological  investigations  for  the  period  July  1,  1953  through 
June  30,  195^«  The  actual  carrying  out  of  the  calculations  indicated 
in  Parts  II  and  III  was  done  under  the  terms  of  the  contract.  The  mathe- 
matical preparations,  i.e.  the  numerical  analysis,  programming  and  coding 
were  carried  out  under  the  terms  of  Contract  No.  N-T-ONE-388,  T.  0.  I 
and  Contract  No.  N-6-OEI-139^  T.  0.  I  between  the  Institute  for  Advanced 
Study  and  the  Office  of  Naval  Research.  Since  the  objectives  of  the 
three  contracts  are  substantially  overlapping  it  was  felt  desirable  to 
include  all  this  material  in  one  report.   In  this  fashion  it  is  hoped 
to  give  the  maximum  possible  information  to  all  interested  agencies. 


John  von  Neumann 
Project  Director 

Institute  for  Advanced  Study 
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PART  I   -  ENGINEEEIWG 


I-l. 


A .   INTRODUCTION 

This  report  describes  the  operation  of  and  engineering  improve- 
ments on  the  electronic  computer  at  the  Institute  for  Advanced  Study 
during  the  period  from  1  July  1953  to  30  June  195^*  The  engineering 
discussion  is  restricted  to  those  features  of  the  present  machine  which 
were  added  during  this  period.  A  full  technical  description  of  the 
machine  prior  to  this  period  is  given  in  the  final  reports  on  Contract 
No.  W-36-03^-OED-7l+8l,  Contract  No.  DA-36-03^-OED-19  (Project  No.  TB3- 
0007F)  and  Contract  No.  DA-36-03^-ORD-1023  (Project  No.  TB3-0007) . 
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B.   SUMMARY 

During  the  first  month  of  the  period  covered  by  this  report  the 
computer  was  in  full  scale  use  at  its  original  construction  location. 
It  was  then  moved  to  its  permanent  location  in  the  mathematical  wing  of 
the  building.  Advantage  was  then  taken  of  the  partial  disassembly  re- 
quired by  the  move  to  carry  out  several  mechanical  and  electrical 
improvements  on  the  memory  unit.   In  addition,  the  control  system,  power 
supplies,  and  cooling  system  were  organized  to  make  operation  and  main- 
tenance much  easier. 

Full  scale  operation  was  resumed  in  December  1953  at  the  new  lo- 
cation.  In  subsequent  months  the  magnetic  drum  was  brought  into  a  state 
of  control,  the  major  work  was  finished  on  a  graphical  output  device  and 
several  other  engineering  improvements  were  made. 


1-3. 


C.  ENGINEERING  WORK 

1.  Machine  move. 

The  computer  was  constructed  in  the  experimental  area 
adjacent  to  the  engineering  laboratory  and  the  shop,  and  after  its 
completion  was  operated  in  this  area  through  all  of  1952  and  the  first 
half  of  1953.  Considerably  operating  experience  was  gained  during  this 
period  and  the  machine  itself  was  brought  into  an  excellent  state  of 
engineering  control. 

However,  the  restricted  space  available  made  operation  cumber- 
some at  best  and  placed  the  machine  and  the  operator  in  the  same  en- 
vironment, this  latter  resulting  in  a  compromise  between  machine  cool- 
ing and  human  comfort.  A  large  space  had  been  provided  in  the  mathe- 
matical wing  of  the  building  and  it  was  accordingly  decided  to  move  the 
machine  into  this  space  and  to  arrange  operating  facilities  as  con- 
veniently as  possible.  An  earnest  attempt  was  made  to  simplify  operat- 
ing procedures,  including  turn  on  and  off,  so  as  to  make  the  machine 
available  to  non- technical  users. 

a.  Memory  improvements  during  move. 

The  overall  height  of  the  machine  was  too  great  to 
permit  passage  through  doorways  so  it  was  necessary  to  separate  the 
memory  and  the  arithmetic  units.  This  occasion  was  used  to  make  several 
improvements  in  the  memory  unit  which  experience  had  shown  desirable: 
i.  Deflection  bus  re -insulation. 

Because  of  the  already  high  average  voltage 
level  of  the  deflection  bus  system  (+13OO  v)  and  the  possibility  of  a 
further  increase  to  improve  read-around,  the  old  fish  paper  insulation 
was  removed  and  nylon  base  phenolic  installed.  In  addition,  teflon 
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insulated  wire  vas  used  for  the  individual  bus-tosocket  connections, 
ii.  Control  plugboard. 

The  focus,  astigmatism,  and  beam  current 
controls  for  each  tube  were  formerly  located  with  the  stage  and  below 
it  with  the  result  that  memory  adjustments  were  executed  from  a  cramped 
position  and  further  involved  stage  by  stage  relocation  of  the  test 
oscilloscope  if  all  stages  were  tuned.   It  was  also  clear  that  the 
potentiometers  themselves  were  causing  day  to  day  variations  in  the 
memory  adjustment,  partly  because  of  their  necessarily  high  resistance 

values.  Therefore  all  cathode  ray  tube  controls  were  moved  to  a  panel 

-39 
at  the  2    end  of  the  machine  and  only  the  individual  electrode  de- 
coupling filters  left  at  the  cathode  ray  tube  sockets.  Further,  the 
focus  and  astigmatism  potentiometer  ensembles  were  abandoned  altogether 
and  replaced  by  a  single  multi-tapped  low  resistance  bleeder  with  a 
plugboard  arrangement.  The  focus  plugboard  supplies  20  voltage  steps 
three  volts  apart,  each  step  feeding  a  row  of  ^1  paralleled  banana 
jacks.  Each  cathode  ray  tube  focus  supply  wire  terminates  in  a  banana 
plug  which  can  be  put  into  any  one  of  the  20  steps.  The  astigmatism 
plugboard  is  similar  except  that  it  provides  12  steps  in  increments  of 
six  volts.  Potentiometers  were  retained  for  the  intensity  control  but 
reduced  in  resistance  from  150  K  to  20  K  ohms.  See  figure  (1)  and  (2). 
iii.  Collected  amplifier  outputs . 

Since  at  least  preliminary  memory  stage 
adjustments  are  made  according  to  observed  dot  and  dash  waveforms,  a 
resistance  divider  was  installed  in  each  discriminator  to  feed  an  at- 
tenuated but  low  impedance  signal  via  a  cable  to  a  kO   position  switch 
at  the  plugboard  location.  Thus  a  complete  preliminary  alignment  can 


Figure  1  -  I-5. 
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be  made  quickly  and  conveniently.  Final  adjustments,  made  according 
to  read -around  test  results,  are  equally  convenient  in  that  a  standard 
focus  or  astigmatism  increment  is  obtained  just  by  shifting  the  plug 
one  row  for  the  stage  in  question.  The  circuit  is  shown  in  figure  (3). 
iv.  Other  changes. 

Additional  changes  were  the  punching  of 
holes  in  the  cathode  ray  tube  support  web  to  facilitate  cooling  air 
flow,  provision  for  a  forty-first  memory  stage  (to  be  completed  if 
needed) ,  consolidation  of  memory  heater  transformer  banks  and  automatic 
screen  charging  at  turn-on. 

b.  AC -DC  turn-on  system. 

A  major  consideration  in  furnishing  the  new  loca- 
tion for  the  machine  was  to  simplify  turning  on  and  turning  off  pro- 
cedures so  that  the  computer  could  be  used  by  the  mathematicians  and 
coders  without  engineering  assistance.  This  would  make  available  the 
i^-S  to  88  hours  per  week  not  regularly  scheduled  and  staffed  by  engineer- 
ing operators.  Accordingly,  a  one  button  system  was  visualized  in  which 
the  single  command   "turn  on"  would  initiate  an  automatic  sequence  of 
events  resulting  in  a  "ready  to  compute"  signal.  Similarly  one  action 
would  be  required  to  "turn  off". 

This  very  simple  system  was  elaborated  slightly  in  the  realiza- 
tion. The  authorization  to  use  the  machine  was  embodied  in  the  require- 
ment of  a  key  to  enable  the  "turn  on"  switch  and  the  following  automatic 
sequence  was  terminated  at  the  completion  of  all  AC  functions.  One  ad- 
ditional button  was  provided  to  apply  DC  voltages  to  the  machine.  A 
"DC  off"  button  was  also  furnished  for  servicing  and  standby,  but  normal 
"turn  off"  was  left  a  one-switch  function. 


Figure  3  -  1-8. 


viDtOA>4P. 
_OOTPUT^ I 


L  J^i^'L'^'!!^'^ 


110 


^      40f>o».Svy/ltc^ 


COLLECTEIO  AMPLIFIEB   OOTPUTS 


1-9. 


Since  the  vacuum  tube  heaters  contribute  about  one -half  of  the 
heat  output  of  the  machine,  it  is  quite  essential  that  the  first  turn- 
on  step  be  to  turn  on  the  cooling  system  and  also  that  continuation  of 
the  turn-on  sequence  be  dependent  on  the  blowers  coming  on.  Following 
this  all  memory  and  arithmetic  unit  heaters  are  turned  up  slowly  by  a 
motor  driven  Powerstat  bank.   In  the  same  period  all  power  supplies  are 
turned  on,  the  large  (PEC)  units  having  their  own  sequencing  arrange- 
ments. Completion  of  these  steps  lights  an  "AC  on"  light.  Turning 
the  DC  on  is  made  a  separate  operation  to  insure  that  someone  is 
present  to  observe  the  machine  at  this  critical  time.  After  IC  turn- 
on  a  time  delay  arrangement  maintains  all  memory  beam  currents  at  a 
higher  than  nornial  value  for  about  two  minutes  in  order  to  bring  the 
charge  on  the  phosphor  to  equilibrium.  When  this  delay  is  over,  the 
machine  is  ready  for  use. 

A  step  by  step  description  of  the  turn-on  process  follows. 
Eeference  is  made  xo  figures  (k)   and  (5) • 

The  starting  control. 

The  whole  machine  can  be  turned  on  only  when  a  key  switch 
is  operated.  This  switch  has  two  functions.  After  the  "AC"  switch  is 
turned  on; 

Turn  key  clockwise  -  energizes  the  air  conditioners  and  blowers. 

Eeturn  key  to  normal  position  -  enables  the  rest  of  the  starting 
process  if  either  compressor  and  booh  blowers  are  properly  energized. 
This  prevents  unnecessary  starting  cycles  in  case  the  air  conditioner 
is  out  of  function  as  the  machine  cannot  be  operated  without  air  con- 
ditioning for  sustained  periods  of  time. 

After  the  key  has  been  returned  to  its  normal  position;  A  limit 


Figure  k   -  I-lo. 
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switch  on  the  Powerstat  driving  motor  senaea  the  poaition  of  the  Power- 
stat  rotor.  If  the  position  of  the  rotor  is  other  than  the  down  limit, 
the  motor  is  energized  so  that  the  down  limit  is  eventually  reached. 
At  that  instant,  due  to  the  operation  of  the  limit  switch,  the  motor 
is  stopped  and  K,  and  Kl^  contactors  are  energized.  K,  is  the  main 
contactor  that  energizes  the  three  phase  AC  supply  to  all  heater  cir- 
cuits. K-,  connects  all  memory,  input-output  organ  and  regulator  heater 
circuits  to  the  output  of  the  Powerstats  and  starts  the  motor  so  that 
the  heater  voltages  are  slowly  raised.  VThen  the  upper  limit  is  reached, 
the  "up"  limit  switch  operates  which  stops  the  motor  and  energizes  K 
which  in  turn  drops  out  K,  and  energizes  K  .  This  act  switches  the 
memory  heaters  from  the  Powerstats  to  Sola  regulating  transformers. 
The  main  (PEC)  power  supplies  are  turned  on  when  K^  is  operated.  The 
PEC  supplies  contain  their  alow-raise  circuits.  At  the  end  of  raise, 
a  relay  operates  to  indicate  completion.  The  release  of  K  reverts 
the  control  of  the  Powerstat  motor  back  to  the  "down"  limit  switch 
which  starts  the  motor  on  the  downward  trip. 

The  Powerstat  is  a  voltage  divider.  The  two  fixed  ends  are 
intrinsically  indistinguishable.  Regardless  of  the  direction  of  travel, 
the  voltage  between  one  of  the  fixed  ends  of  the  divider  and  the  rotor 
is  increasing.  This  fact  is  used  to  raise  the  heater  voltages  of  the 
arithmetic  units  by  the  action  of  K-,. 

At  the  end  of  the  downward  travel,  the  "down"  limit  switch 
stops  the  motor  and  transfers  from  K^  to  &  whereby  the  arithmetic 
heaters  are  now  supplied  directly  from  the  three  phase  AC  lines. 

Notice  that  after  K  has  been  energized,  the  air  conditioning 
monitoring  circuit  no  longer  has  control  of  the  machine.  This  is  to 
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prevent  interruption  of  machine  operation  by  normal  cycling  of  the  com- 
pressors (under  thermostat  control) . 

The  DC  turn-on  is  conditional  on  the  PEC  supplies.   If  these  sup- 
plies are  in  operating  condition,  the  DC  turn-on  contactors  can  be  enabled. 
A  manual  PEC  disable  switch  is  installed  to  facillitate  servicing. 

In  figure  (k)   the  AC  starting  circuitry  are  located  on  the  left  of 
terminal  9>  while  on  its  right  are  DC  switching  circuitry.  The  DC  switch- 
ing are  arranged  such  that  the  input-output  organ  and  drum  units  can  be 
separately  turned  off  while  the  machine  DC  is  on,  to  enable  maintenance 
on  those  units . 

In  all  cases  of  on-off  procedure,  the  -300-volt  bus  is  switched 
on  one  relay  time  later  than  all  plus  voltages.  This  is  to  prevent 
excessive  heater -cathode  voltages  in  the  absence  of  plate  voltages  on 
cathode  followers  and  coincidence  circuits.   On  turn-off  the  reverse 
is  made  true.  The  heaters  are  also  prevented  from  being  turned  off 
while  DC  is  applied  to  the  machine. 
c.  The  new  location, 
i.   Layout. 

Figure  (6)  shows  the  physical  layout  of  the 
new  room.  The  computer  is  in  a  separate  room  from  the  operator,  but  can 
be  viewed  through  a  large  window  in  front  of  the  operating  table.   The 
computer  room  itself  need  not  be  entered  except  for  maintenance  and 
trouble  shooting.   It  can  be  kept  cooler  than  human  comfort  temperature 
and  also  more  free  of  dust  than  if  operators  were  in  the  same  room. 

The  operating  room  contains  all  the  needed  operating  controls, 
a  monitor  cathode  ray  tube  for  displaying  the  contents  of  any  memory 
stage,  the  graphing  device  (described  elsewhere  in  this  report),  and 
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an  IBM  type   51^  reproducing  punch  which  is  used  for  input  and  output. 
A  third  room  houses  the  IBM  equipment  needed  for  producing  the 
original  machine  input  decks  and  for  processing  outputs.   It  contains  two 
alphameric  key  punches,  a  reproducing  punch,  an  alphameric  tabulator  and 
a  sorter. 

ii.  Power. 

To  meet  the  power  requirements  of  the  computer 
and  its  associated  equipment  a  200  ampere  feed  was  installed  from  the  main 
building  load  center  to  the  machine  location.  Here  the  load  was  divided 
into  three  parts ; 

1)  the  high  current  DC  supplies  -  6o  amperes  per  phase, 

2)  blowers  and  cooling  system  -  ko   amperes, 

3)  machine  vacuum  tube  heaters  and  miscellaneous  -  50  amperes. 
The  switching  of  these  loads  during  normal  machine  turn-on  and  off  is 
controlled  by  the  starting  control  (see  C  1.  b.).  There  are  four  of  the 
high  current  DC  (PEC)  supplies  which  carry  the  major  DC  machine  load 
directly,  but  some  power  is  also  needed  at  intermediate  voltage  levels. 
These  are  provided  by  a  bank  of  secondary  regulators  located  in  the  power 
switching  center.   (See  C  2.  b.)  All  power  runs  to  the  machine  are  in 
conduits  under  the  floor  or  in  the  ceiling.  As  a  precaution,  wire  screen 
was  laid  on  the  floor  and  ceiling  and  both  interior  walls  to  form  an 
open  sided  box  enclosing  the  machine.  All  conduits  were  bonded  to  this 
screen  at  the  point  of  entry  and  the  machine  base  bonded  to  floor  screen. 

The  starting  control,  regulating  transformers,  secondary  regu- 
lators, and  the  high  voltage  supplies  for  the  memory  are  all  mounted  in 
a  steel  frame  and  constitute  the  power  switching  center.  This  is  shown 
in  figures  (7)  and  (8) . 
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ill.  Cooling. 

A  closed  circuit  air  cooling  system  was  in- 
stalled in  order  to  provide  clean,  low  hximidity  cooling  air  to  the  ma- 
chine. Air  is  blown  through  a  floor  duct  into  the  base  of  the  computer, 
rises  through  it  and  exhausts  through  a  ceiling  duct,  returns  through 
an  exhaust  blower  air  filters  and  cooling  coils  to  the  floor  duct  again. 
Cooling  is  provided  by  two  remotely  located  7  l/2  ton  compressors  and 
h.eat  is  dissipated  by  a  Unicon  refigerant  to  air  heat  exchanger  located 
on  the  building  roof.  The  refrigerant  system  of  each  compressor  is  sep- 
arate so  that  maintenance  can  be  done  on  one  system  while  the  other  tem- 
porarily carries  the  machine.  Because  this  cooling  system  is  intended 
for  year-round  operation,  pressure  regulating  bypass  valves  were  installed 
on  the  Unicon  to  maintain  constant  compressor  head  pressures  regardless 
of  outside  air  temperature. 

The  air  circulation  path  is  shown  in  figure  (9)  and  the  power 
wiring  in  figure  (10). 
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2.  Engineering  Improvements . 

A  number  of  changes  and  improvements  have  been  made  on  the 
machine  to  increase  ease  and  versatility  of  operation  and  also  reliability. 
These  changes  are  set  forth  in  this  section  as  a  series  of  articles, 
a.  RI  digit  insertion, 
i.  General. 

A  considerable  operating  convenience  on  the 
machine  is  a  means  for  manually  inserting  words  in  the  memory,  as,  for 
instance,  parameter  changes  between  one  run  and  the  next.  Since  storage 
into  the  memory  is  from  RI  (the  accumulator  register) ,  it  suffices  to 
provide  a  means  for  modifying  the  contents  of  RI  by  a  switch  bank  at  the 
operating  position. 

Several  methods  of  remote  toggle  switching  have  been  devised;  some 
of  them  have  been  used  at  this  project.  One  method  is  to  pull  on  the 
proper  plate  of  a  toggle  by  means  of  a  voltage  impressed  on  a  wire  to 
that  plate.  A  change  is  enforced  in  the  level  of  voltage  seen  at  the 
plate,  and  the  toggle  will  flip.  This  process  is  wasteful  of  power,  re- 
quires local  decoupling  for  optimum  noise  reduction,  imposes  a  large 
constant  capacitative  load  on  the  toggle  network,  and  is  in  general  in- 
elegant. Further,  it  is  not  simple  to  get  bidirectional  control  with 
only  one  wire  in  a  plate  control  scheme. 

If  a  grid  could  be  easily  swung  up  and  down  by  a  controlling 
voltage,  we  should  be  dealing  with  information  control  rather  than  power 
control,  and  one  would  simply  have  to  apply  one  of  two  voltages  on  a  line 
to  obtain  either  toggle  state  at  will.  Obviously  it  is  not  desirable  to 
burden  a  grid  control  point  with  an  appreciable  load.  Loss  of  operating 
speed  and  reliability  would  be  an  inevitable  consequence.   It  is  possible 
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to  uae  diode  or  triode  elements  placed  close  to  the  controlled  point  to 
act  as  decoupling  agents  for  the  long  controlling  line.  This  sort  of 
solution  is  neither  economical  nor  elegant. 

A  small  neon  lamp,  such  as  the  HE-2  has  interesting  properties 
which  lend  themselves  to  this  decoupling  problem.  In  the  off  state,  a 
typical  lamp  of  this  type  has  a  capacity  between  its  leads  of  t^   O.7  uuf 
and  a  resistance  on  the  order  of  several  thousand  megohms.  After  being 
turned  on,  a  typical  NE-2  has  a  resistance  on  the  order  of  65  K  with  one 
ma  through  it  and  35  K  when  passing  2  ma. 

Figure  (11)  shows  a  possible  use  of  the  NE-2  for  toggle  control. 
During  the  time  that  no  control  is  being  exercised,  the  point  x  on  the 
toggle  has  incurred  no  more  loading  than  a  capacity  <  O.7  uuf.  This  is, 
for  all  practical  purposes,  a  trivial  perturbation.  The  bias  existing 
on  the  control  line  at  this  time  is  chosen  to  be  a  value  halfway  between 
the  ordinary  extremes  of  toggle  voltage  at  point  x.  This  insures  thiat 
neon  firing  will  not  occur  accidentally  (since  x  commonly  has  a  stretch 
of  a»  ^0  v.) .  Most  NE-2  samples  will  not  fire  at  less  than  75  applied 
volts,  so  we  are  doubly  safe. 

Switching  either  to  E  or  E  then  pulls  x,  the  right  toggle  grid 
to  a  point  where  stability  oecijirsj  on  returning  the  switch  to  the  neutral 
or  bias  position,  the  new  toggle  state  remains. 

For  the  condition  of  switching  to  E  ,  grid  ciirrent  will  flow  in 
the  right  side  of  the  toggle  due  to  the  grounded  cathode.  The  values 
of  E  and  E  are  chosen  to  limit  this  current  to  a  safe  value.  A  range 
of  from  1.5  to  1.75  ina  has  been  found  to  provide  certain  yet  safe  opera- 
tion for  these  toggle-neon  combinations. 

A  relatively  simple  extension  of  circuitry  then  allows  flexible 
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control  of  a  number  of  such  stages,  single  or  multiple  switching  then 
being  possible,  as  seen  in  figure  (12). 

S,  ia  seen  to  be  the  individual  stage  control  while  S  will  simul- 
taneously switch  any  number  of  stages  attached  to  the  common  bus.   It  is 
desirable  to  apply  E  bias  in  such  a  way  as  to  insure  its  application  at 
all  times,  not  depending  on  a  set  of  closed  switch  contacts.   Somewhat 
more  provision  against  undesirable  noise  voltages  appearing  on  the  con- 
trol lines  is  obtained  this  way. 

Modifications  of  the  above  system  are  in  use  at  this  project  and 
have  been  found  to  operate  satisfactorily.  A  simple  arrangement  of  com- 
mon components  provides  remote  toggle  control  with  great  economy  of 
space,  power,  and  cost. 

ii.  Detailed  description. 

A  test  of  ten  KE-2  miniature  neon  lamps 
selected  at  random  from  new  stock  was  conducted  to  obtain  conduction, 
firing  and  extinction  characteristics.  The  results  are  tabulated  on 
the  following  page. 

Of  the  ten  lamps  tested  the  highest  voltage  required  for  ignition 
was  86  V  and  the  lowest  voltage  which  would  insure  sustained  conduction 
before  extinction  was  57  v.  Thus  for  an  unselected  group  of  IIE-2's,  we 
must  insure  >  86  volts  for  firing,  <  57  volts  for  extinction,  and  for 
nominal  currents  through  the  lamp  during  conduction  the  voltage  must  be 
between  ~  68  and  7^  volts. 

A  graphic  suimnary  of  these  characteristics  is  seen  in  figure  (13). 

In  October  1953  a  switching  facility  incorporating  these  lamps 
was  built  and  installed  as  a  control  for  the  E,  register  of  the  machine. 
A  bank  of  8l  switches  and  associated  resistors  was  provided  in  a  control 
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M.A.  THROUGH  LAMP 
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box  located  at  the  machine  operating  position.  There  are  ko   switches 
for  putting  ones  and  kO   for  putting  zeros  into  the  register.  Another 
switch  allows  the  entire  register  to  he  cleared  either  to  one  or  to 
zero.  The  circuit  for  this  control  unit  is  shown  in  figure  (1^+). 

A  standby  bias  potential  of  22-5  v  is  applied  to  the  side  of 
the  neons  remote  from  the  controlled  toggle.  Thus  at  no  time  will  a 
voltage  differential  of  >  25  v  appear  across  the  NE-2  and  accidental 
firing  is  inhibited. 

The  33  K  resistors  in  the  switching  lines  provide  sufficient  cur- 
rent limiting  to  protect  both  toggles  and  neons.  In  the  event  of  the 
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worst  combination  of  !(%   toggle  resistor  drifts,  less  than  3  ma  of  grid 
ciirrent  will  flow  when  the  toggle  is  cleared  to  the  one  state.  The  sup- 
ply voltages  of  +  135  volts  are  sufficient  to  provide  safe  operation  for 
any  combination  of  10^  resistor  drifts. 

This  system  has  been  in  use  for  over  six  months  and  has  given 
satisfactory,  reliable  operation.  The  neon  lamp  used  as  a  decoupling 
control  device  seems  to  be  a  feasible  solution  to  a  large  class  of  com- 
puter operation  problems,  and  will  undoubtedly  see  more  service  in  the 
present  machine. 

b.  New  secondary  regulators. 

The  principal  DC  load  of  the  machine  flows  from  the 
high  current  +  110  volt  and  +  2ii-0  volt  supplies  and  returns  through  either 
ground  or  the  -300  volt  supply.  However,  smaller  amounts  of  power  are 
needed  at  other  voltage  levels,  some  of  them  higher  than  +  2*4-0  so  a 
fourth  high  current  supply  is  at  +  380  volts  to  provide  a  memory  volt- 
age directly  and  also  to  feed  a  number  of  secondary  regulators  which  in 
turn  produce  the  needed  intermediate  levels,  particularly  to  the  adder 
and  digit  resolver. 

A  set  of  such  secondary  regulators  had  been  constructed  along 
with  the  arithmetic  unit  and  were  continued  in  use  through  the  machine 
move.  These  were  essentially  independent  units,  each  with  its  own 
heater  supply,  switching,  metering,  and  voltage  reference,  mounted  in 
relay  racks  with  the  tubes  horizontal.  They  were  unsatisfactory  for 
several  reasons:  the  two  stage  error  amplifier  used  was  only  condi- 
tionally stable;  the  independent  construction  led  to  excess  bulk  and 
failure  hazard;  and  the  horizontal  mounting  of  the  regulator  tubes 
gave  considerable  grid-cathode  short  trouble.  Therefore  a  new  regulator 


1-28. 


system  was  designed  and  constructed  to  meet  the  now  definite  power  needs 
of  the  machine. 

The  new  system  consists  of  two  15"  x  33"  hinged  horizontal  chassis 
side  by  side  with  all  tubes  mounted  vertically.  Each  chassis  accommo- 
dates only  tubes  and  resistors,  with  all  filament  transformers  and  filter 
condensers  mounted  on  a  shelf  below.  Fewer  tubes  are  required  in  the  new 
system  because  of  shunts  across  the  6o80's.  These  shunts  supply  the 
steady  component  of  current.  A  single  reference  chassis  is  provided  for 
all  regulators  (see  figure  (15)) •  The  new  regulators  were  installed  in 
May  195^  and  to  date  have  proven  to  be  quite  trouble-free  in  that  there 
has  been  only  one  failure  --a  microphonic  6o8o.  The  installation  was 
accompanied  by  a  reduction  in  the  number  of  voltage  levels  required  for 
machine  operation  and  a  program  for  further  reduction  is  under  study. 

A  typical  regulator  is  shown  in  figure  (l6) .  The  elements  include 
a  6AU6  pentode -amplifier  and  one  or  more  6o8o's  as  the  series  tube(s) 
dependent  on  load  requirements.  The  6o80's  are  operated  at  a  maximum  of 
half  of  rated  plate  dissipation  and  half  rated  zero  bias  emission,  deliver- 
ing a  maximum  of  about  110  ma.  at  55  volts.   In  series  with  each  section 
of  every  6o8o  are  two  10-ohm  resistors,  the  one  nearest  the  plate  serving 
as  a  suppressor  and  the  other  providing  a  means  of  observing  the  current 
flow  in  that  section  (suitable  wiring  and  tip- jacks  allowing  these  voltage- 
drop  readings  to  be  taken  remotely) .  By  periodic  observation  of  the  pro- 
portionality between  the  readings  for  all  paralleled  sections  we  can 
evaluate  the  relative  performance  of  the  series  tubes  and  make  any  re- 
placements as  required.  As  yet  it  has  not  been  necessary  to  make  any 
replacements  except  for  the  6o8o  mentioned  above. 

The  output  voltage  of  each  regulator  can  be  varied  about  10  volts 
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above  or  below  its  nominal  value  and  will  regulate  to  within  no  more 
than  one  volt  for  all  normal  loads  over  thia  range.  DC  voltage  inputs 
for  all  of  the  regulators  and  for  the  reference  chassis  are  derived  from 
commercial  high-current  regulated  power  supplys  which  also  feed  the  com- 
puter directly. 

The  various  voltage  levels  and  protective  devices  associated  with 
the  machine  are  not  automatically  monitored  at  present  but  such  a  system 
is  in  the  initial  stage  of  development.  This  system  will  detect  and 
record  momentary  or  fixed  variations  in  all  DC  levels  and  also  indicate 
the  state  of  the  input  fuses  to  the  various  regulators  and  power  supplies. 
Any  blown  fuse  will  cause  all  DC  to  be  removed  from  the  machine,  but  for 
any  DC-level  variation  which  exceeds  the  normal  there  will  merely  be  an 
alarm.  The  indications  of  all  failures  so  observed  will  remain  avail- 
able in  the  absence  of  all  machine  voltages. 
c-      Present  drum  techniques . 

The  drum  which  has  been  described  in  previous 
reports  is  now  used  extensively  by  practically  all  coders.  Certain 
engineering  changes  have  been  incorporated  into  its  circuitry  with  the 
intention  of  obtaining  more  trouble-free  operation.  Among  these  are 
a  reduction  in  the  number  of  relays  and  tubes  in  the  control  and  also 
relocation  of  some  elements  so  as  to  facilitate  preventive  'flaaintenance. 
The  surface  coating  is  now  erased  to  magnetic  neutral  rather  than  to  a 
polarized  state  as  was  previously  described.  A  square-wave  generator 
is  used  at  present,  the  amplitude  being  reduced  slowly  by  manual  ad- 
justment, but  shortly  there  will  be  facility  for  erasing  automatically 
simply  by  pushing  an  "ERASE"  button.  Most  of  this  circuitry  now  exists 
and  when  operative  will  function  as  follows:  A  unique  marker  (Sync  1-2) 
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will  trigger  a  one-shot  multivibrator  and  for  one  complete  revolution 
the  drum  will  be  erased  "negatively"  until  the  arrival  of  Sync  1-2  again 
at  which  time  the  multivibrator  will  provide  for  "positive"  erasing. 
This  alternation  persists  for  about  90  seconds  during  which  time  the 
erasing  current  is  continually  reduced  by  means  of  an  E-C  time  constant 
imposed  on  a  6o8o  tube  in  series  with  all  ko   pulsers.  Erasing  is  now 
done  periodically  (at  least  once  a  week)  and  it  appears  that  this  pro- 
cedure reduces  the  occurrence  of  errors, 
d.  Graphing  circuit. 

The  design  of  the  graphing  function  for  the  machine 
was  initiated  early  in  April  195^-  This  function  is  intended  to  provide 
an  auxiliary  facility  for  maintaining  a  check  on  the  progress  of  a  com- 
putation.  It  is  expected  that  considerable  time  and  effort  can  be 
saved  by  having  available  a  means  for  visually  observing  some  of  the 
variables  of  a  solution  while  the  computation  is  actually  being  performed. 
One  may  see  continuity  or  discontinuity  in  a  function,  for  instance,  or 
observe  the  number  of  oscillations  between  two  given  points.   It  is  use- 
ful also  to  know  whether  two  processes  being  initiated  from  divergent 
points  tend  to  approach  each  other  and  meet  in  a  smooth  fashion.  A 
feeling  for  scaling  quantities  as  well  may  be  obtained  from  such  a  dis- 
play. These  various  factors  are  obtainable  without  such  a  visual  aid 
only  by  means  of  tedious  manual  collation  and  computation.  Of  further 
use  would  be  the  photographing  of  graphed  displays  for  future  analysis. 
The  display  is  designed  to  be  observed  on  a  7"  cathode  ray  tube. 
This  tube,  a  TJPl,  has  a  medium  persistence  green  trace  and  is  similar 
to  the  memory  and  slave  tubes  used  at  present  in  the  main  machine. 
Facility  for  switching  either  the  graphing  or  the  slave  function  to  this 
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tube  is  provided,  so  that  when  graphing  is  not  desired  an  auxiliary, 
enlarged  memory  slave  monitor  is  available.  Diiring  such  slave  monitor- 
ing (when  it  is  used  in  conjunction  with  the  existing  memory  slave  unit) 
it  becomes  possible  to  observe  two  separate  memory  stages  at  once,  or  to 
monitor  two  pairs  of  "compared"  stages.  The  resolution  of  the  graphing 
display  is  9  binary  digits  or  512  points  in  both  horizcntal  and  vertical 
coordinates.  Reasons  for  this  choice  will  be  given  when  the  details  of 
the  deflection  system  are  described.  Kesolution  during  memory  slave 
operation  is,  of  course,  the  standard  raster  of  32  x  32  spots. 

Information  input  to  the  display  is  obtained  from  two  sources. 
During  graphing  we  have  access  to  the  output  of  the  forty  digits  of  the 
magnetic  drum  (Group  "0"  or  "1"  heads) .  We  wish  to  use  two  groups  of 
nine  digits  each,  one  for  horizontal  and  one  for  vertical  deflection. 
By  suitable  selection  from  this  kO   digit  input  it  is  possible  to  ob- 
tain digits  of  properly  selected  weight  and  position.  For  slave  opera- 
tion, the  digit  input  comes  directly  from  the  ten  digit  address  of  the 
machine's  dispatch  counter  (center  rank). 

i.  Deflection  for  graphing  operation. 

Forty  digits  are  available  from  the  output 
of  the  magnetic  drum  from  either  of  the  two  groups  of  heads.  Either 
group  (zero  or  one)  is  selected  at  the  will  of  the  operator  and  the  forty 
digits  are  then  fed  from  the  drum  over  a  cable  to  kO   corresponding  triodes 
in  the  graphing  unit.  These  digit  lines  have  origins-ted  from  the  drum 
head  amplifiers  which,  with  the  IBM  output  circuitry  are  wired  to  cathodes 
of  triodes  at  the  main  machine's  input  for  these  signals.   It  is  necessary 
that  each  of  these  digit  lines  be  isolated  from  the  machine  during  ordinary 
computation  and  that  they  be  connected  for  an  external  order.  The 
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requirement  then  is  that  for  computing  theae  lines  be  at  >+  5  v.  and 
during  an  external  order  they  fall  to  >  -  25  v.  A  condition  which  is 
unique  at  these  times  is  the  presence  or  absence  of  lEG  ,   and  this  signal 
may  be  used  to  change  the  state  of  these  lines.  The  grids  of  each  of  the 
14-0  digit  input  line  triodes  on  the  graphing  chassis  therefore  are  re- 
turned to  a  common  bus  which  is  supplied  by  a  driver  system  under  the 
control  of  ING^-  During  an  external  order  then  these  grids,  hence  their 
cathodes  are  down  and  the  digit  lines  are  in  communication  with  the  com- 
puter. At  all  other  times  the  grids  being  up,  pull  their  respective 
cathodes  to  «<  +  5  v.  which  insures  the  isolation  of  the  lines  from  the 
machine  yet  retains  the  function  of  transfer  of  information  from  the 
drum  amplifier  outputs  to  graphing  input. 

The  outputs  of  the  digit  line  triodes  are  wired  to  a  circular 
array  of  connector  jacks.  Two  "combs"  of  plugs  may  then  be  inserted 
into  these  jacks  so  that  two  sets  of  nine  consecutive  digits  each  may 
be  picked  off.  Thus  it  is  possible  to  select  the  I8  digits  to  be  used 
from  among  the  ii-O  available  in  order  to  obtain  the  desired  position  and 
"weighting"  of  digits.  An  alternative  patch  plug  arrangement  is  pro- 
vided so  that  one  may  choose  other  than  consecutive  sequences  of  digits. 

The  selected  signals  (nine  for  horizontal  deflection  and  nine 
for  vertical)  are  then  applied  through  relay  contacts  (the  upper  level 
of  these  signals  being  diode  bumped  at  ground)  into  gates.  At  the  proper 
moment  in  time  this  information  is  gated  into  a  set  of  toggles  which  feed 
the  deflection  circuitry. 

Each  toggle  (one  for  each  of  the  18  digita)  controls  the  current 
in  one  or  the  other  side  of  a  dual  triode  deflection  adder  tube  which 
supplies  two  deflection  summing  buses.  This  is  essentially  the  scheme 
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uaed  in  the  computer  Williams  memory  system  anfl  is  described  in  detail 
in  a  previous  report. 

There  are  now  two  pairs  of  summed,  balanced  deflection  corrents, 
one  pair  for  horizontal  and  one  for  vertical,  which  are  passed  through 
the  deflection  driver  tubes.  Voltages  developed  across  the  plate  re- 
sistors of  these  tubes  are  then  directly  applied  to  the  deflection 
plates  of  the  cathode  ray  tube. 

The  choice  of  a  512  x  512  array  was  arrived  at  from  the  following 
considerations.  The  usable  projection  area  of  a  7"  cathode  ray  tube  is 
roughly  that  defined  by  a  6  inch  diagonal  or  a  square  with  k  l/k   inch 
sides.  One  thousand  and  twenty-four  spots  along  a  k  l/k   inch  line  would 
have  a  spacing  of  a;  4  mils  while  51^  spots  would  be  spaced  ~  8  mils 
apart.  Considering  the  currents  necessary  in  the  adder  tubes,  let  us 
assume  that  25  m.a.  is  a  reasonable  maximum  current  for  the  largest  de- 
flection digit.  Then  for  1024  spots  obtained  from  10  binary  digits,  the 
smallest  deflection  current  would  be  -^  ' — -   (where  n  =  10)  or  approxi- 
mately 50  u.a.  For  512  spots,  the  minimjAm  digits  current  ia  100  u.a. 
The  amoxmt  of  noise  and  uncertainty  which  would  be  fouid  in  any  practical 
circuit  is  not  inconsiderable  compared  to  50  u.a.,  and  requiring  a  mini- 
ature tube  to  pass  as  much  as  50  m.a.  is  not  a  simple  task.  Since  the 
e3re  cannot  easily  resolve  500  lines  or  dots  in  a  U  or  5  inch  length  when 
viewed  a  foot  or  two  distant,  it  ia  considered  reasonable  to  design  for 
a  nine  digit  or  512  spot  resolution.  For  all  practical  purposes  a  func- 
tion plotted  with  this  resolution  will  appear  as  a  continuous  line,  and 
with  care  a  circuit  can  be  made  to  function  reliably  with  adder  currents 
ranging  between  25  m.a.  and  100  u.a. 
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ii.  Deflection  for  slave  operation. 

The  five  digits  each  for  horizontal  and 
vertical  deflection  are  brought  over  to  the  graphing  unit  from  the  ma- 
chine's dispatch  counter.  After  passing  through  amplifier  tubes  the 
digits  are  applied  through  relay  contacts  to  the  moat  significant  digit 
lines  entering  the  set  of  deflection  toggles.  Thus  with  the  relays  E^ 
and  E  in  the  slave  position,  the  normal  machine  raster  deflection  is 
produced  by  the  adder  in  the  same  manner  that  the  gi-aphing  raster,  just 
described,  is  produced.  A  32  x  32  spot  pattern  is  then  seen  on  the 
cathode  ray  tube. 

iii.  Beam  turnon. 

Beam  turnon  control  for  the  graphing  opera- 
tion is  derived  from  the  graphing  control  circuit  and  will  be  described 
below  in  the  section  dealing  with  the  control. 

For  slave  operation  the  beam  turnon  signals  are  generated  by  a 
circuit  comprising  a  selector,  comparator,  and  a  stretching  circuit.  The 
ko   signals  are  brou^it  over  cs-bies  from  the  machine's  beam  turnon  circuits 
to  two  ko   position  selector  switches  on  the  graphing  chassis.  The  signals 
thus  selected  (A  and  B)  are  amplified  and  bumped  between  gi'ound  and  -  10. 
Ground  level  represents  beam  turnon.  Each  signal  is  then  introduced  into 
the  channel  comparator-selector  circuit  shown  in  figure  (l8). 

If  the  selector  switch  is  in  the  "A"  position,  the  plate  of  Y^    is 
connected  to  the  output  line.  Vp  is  cut  off  at  this  time  by  the  trans- 
pose network.  The  grid  of  V,  is  pegged  at  «;  -  5  v.   Its  cathode  being 
bumped  at  ground  by  the  diode  on  the  signal  line,  cannot  rise  above  0  v. 
V,  is  then  essentially  cut  off  until  a  beam  turnon  signal  of  -  10  comes 
along.  At  this  time  V^  will  conduct  and  the  pulse  will  be  seen  as  a  down- 
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going  signal  on  the  output  line.  Similarly  channel  "B"  will  conduct 
turnon  signals  with  the  selector  switch  in  the  "B"  position. 

For  comparison  the  outputs  of  V,  and  V  are  connected  together 
and  their  plates  observed.  Each  grid  will  he  at  «%»  -  5  v.  when  each 
cathode  normally  rests  at  ground,  thus  both  tubes  are  cut  off.   If  sim- 
ilar signals  are  applied  to  each  input,  both  grids  and  cathodes  will  dip 
similarly  and  no  output  signal  will  be  developed.  However,  if  dissim- 
ilar signals  are  present,  an  output  signal,  showing  that  a  difference 
exists,  is  developed.  Suppose  signal  "A"  is  present  and  "B"  is  absent. 
The  grid  of  V,  is  at  -  5  v.  while  the  cathode  of  V^  goes  to  -  10.  Y^ 
conducts  and  a  signal  appears  on  its  plate.  The  cathode  of  V^  remains 
pegged  at  ground  but  its  grid  drops  to  -  15  v.  and  thus  V  is  cut  off. 

The  beam  turnon  signals  appearing  on  the  output  line  are  then 
ready  for  transmission  to  the  cathode  ray  tube.  For  slave  operation  the 
relatively  short  times  dxiring  which  the  beam  turnon  signals  exist  (1.2 
us  .  for  normal  dot  and  k.l   us.  for  normal  dash)  are  insufficient  to 
guarantee  maximum  slave  brightness.   If  these  turnon  pulses  could  be 
lengthened,  the  slave  intensity  would  correspondingly  improve. 

A  "boot-strap"  circuit  is  used  to  stretch  the  beam  turnon  pulse. 
The  present,  normal  signals  are  shown  in  figure  (19)' 

The  turnon  pulses  can  be  lengthened  by  a  factor  of  ~  3  to  1  for 
slave    without  exceeding  the  clock  period  limitations. 

The  circuit  used  to  lengthen  the  signals  is  shown  in  figure  (20). 

During  beam  off  time  E.   is  at  +  10  and  Ec  is  bumped  at  +  110, 
its  lower  level  by  grid  current  in  V  .  The  output  at  the  plate  of  V 
at  this  time  is  "down"  or  the  beam  turnon  signal  is  in  the  off  condition. 

When  E  in.  goes  negative  the  left  half  of  V^  is  cut  off  and  Ec 
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rises  from  its  lower  level  of  +110  at  a  rate  determined  by  the  current 

through  R  which  is  charging  C  This  current  I^^  =  (380  -  110) /82  K  =  3.3  ma. 

which  will  charge  a  150  p  f.  capacitor  at  tlrie  rate  of  22  v. /usee. 

When  E  in.  goes  positive,  the  "constant"  cathode  current  of  the 

switching  tube  V,  is  impressed  on  capacitor  C  to  bring  Ec  back  down  to 

its  standby  condition  of  +  llOv.  Depending  on  the  setting  of  potentiometer 

E^  this  cathode  current  will  be; 

7  ma  if  E  =  0 

k.2   ma  if  E^  =  30  K. 

When  Ec  moves  more  than  a  few  volts  above  +  110,  Y^   cuts  off  and 

its  plate  signals  the  CEO  grid  to  turn  the  beam  on.  The  beam  will  remain 

on  until  Ec  returns  close  enough  to  +  110  to  cause  V  to  conduct  again. 

For  our  purposes  here  we  may  assume  that  ¥_  cuts  off  as  soon  as  Ec  begins 

its  rise. 

The  actual  stretching  may  be  considered  as  follows  (see  figure  (21)): 

During  time  T,  :  Ic  =  Ic^^  =  Ij^  ^-^  ^  '^'^ 

Dijxing  time  T^  :  Ic  =  Ic^  =  I^  -  Ig  (note  that 

It   >   Ir   ) 

"^1  ■*"  ^2       '^l/'^2  ^  '■ 
Stretching  ratio 


T^  T^/r, 


but  T^/Tg  =  in       =  It  -% 
I™ 


h.2  -  3.3  7  -  3.3 


J,  2                                7 
the  range  of  = =  i+.7  to  =  1.9 
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iv.  Control. 

The  graphing  operation  ia  very  similar  to  the 
normal  magnetic  drum  operation  in  that  blocks  of  32  words  each  are  selected 
for  use  by  a  ten  digit  "priming  order".  The  first  five  digits  of  this 
order  specify  the  number  of  the  starting  block  while  the  second  five  digits 
(the  right  half  of  the  order)  specify  the  number  of  blocks  to  be  read  out. 
During  normal  drum  operation  this  priming  order  is  derived  from  the  computer 
while  for  graphing  the  ten  priming  digits  are  set  up  by  a  bank  of  switches 
on  the  graphing  control  panel.  As  long  as  a  graphing  cycle  is  in  progress 
(the  machine  not  executing  an  external  order) ,  the  graphing  priming  order 
takes  precedence  over  the  normal  priming  number  in  the  drum  register. 

For  graphing  the  priming  order  which  has  been  set  up  on  the  control 
panel  is  waiting  to  be  transferred  to  and  recognized  by  the  drum's  coinci* 
dence  circuit  input  selector.  The  ten  digits  of  the  order  are  applied  to 
the  coincidence  circuit  input  through  NE-2  neon  lamps  acting  as  isolating 
switches'  The  use  of  these  neons  as  linkage  devices  has  been  described 
in  another  section  of  this  report.  Sufficient  compliance  is  inserted  be- 
tween the  toggles  of  the  priming  register  and  the  coincidence  gates  to 
insure  that  the  remote  insertion  of  priming  digits  does  not  affect  the 
state  of  the  priming  toggles.  This  provision  allows  the  last  regularly 
stored  priming  number  to  remain  intact  for  use  or  reference  at  any  time. 

A  typical  stage  is  shown  in  figure  (22).  If  the  digit  being  in- 
serted into  the  right  half  of  the  gate,  as  shown  is  one  of  the  first  five 
digits  of  the  priming  number,  then  the  left  half  of  the  gate  will  be  as- 
signed to  a  digit  in  the  second  half  of  the  priming  number.  If  the 
second  half  digit  is  being  switched  to  the  negative  switching  potential, 
the  left  half  of  the  gate  will  not  be  affected  during  the  period  of  no 
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gating;  however,  should  the  digit  be  switched  to  the  positive  potential, 
the  absolute  voltage  appearing  at  the  left  gate  grid  will  not  have  a 
guaranteed  value  (due  to  variations  from  one  KE-2  to  another)  and  it  is 
conceivable  that  gating  may  inadvertantly  occur  prior  to  the  desired  gat- 
ing time.  Thus  it  is  necessary  to  delay  the  application  of  the  positive 
potential  at  the  graphing  priming  switches  until  the  time  for  the  use  of 
that  particular  group  arrives.  This  ia  done  by  utilizing  the  BO/3  toggle 
of  the  drum  to  apply  either  of  the  two  positive  switching  voltages  by 
means  of  Vp,  and  '^       shown  in  the  accompanying  control  circuit  schematic, 
figure  (25) . 

To  begin  the  graphing  cycle,  let  us  assume  that  T-^  has  just  been 
turned  on.   It  is  seen  that  the  two  main  control  lines,  line  1  and  line  2 
are  pulled  up  and  down  respectively  by  means  of  triodes  V--  and  V,  .  Con- 
aider  the  functions  controlled  by  line  2:  diode  V  will  pull  down  on  the 
sync  permission  line,  enabling  the  magnetic  drum  sync  signals  to  be  ac- 
tivated. Diode  V_  pulling  down  on  the  T^_,  control  line,  conditions  the 
^omm  toggle  (in  the  magnetic  drum  turn  on  circuits)  so  that  the  very  next 
main  timing  pulse  (sjnic  1,  2)  turns  on  T„,_  initiating  a  graphing  cycle. 
The  seven  digit  counter  in  the  drum  coimts  sync  1  pulses  and  when  this 
count  equals  the  "starting  block  number"  of  the  priming  register  a  counter 
satisfy  signal  turns  on  the  T^^.^^^  toggle.  Meanwhile  (via  line  1)  the 
normal  consequences  of  T-T-nn-c'  being  ON  in  the  drum  control  are  inhibited. 
It  should  be  noted  that  the  coincidence  circuit  has  been  referring  to  the 
left  half  of  the  priming  register  for  this  operation  since  both  grids  of 
V   are  cut  off  —  i.e.  line  2  of  the  control  is  negative  and  '^■dqH   is 
in  the  OFF  state.  This  supplies  a  positive  neon  switching  voltage.  For 
any  digits  requiring  a  negative  neon  switching  voltage  the  following 
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occurs:  V.^o  is  conducting  during  this  tirne  s-;ince  line  1  is  up,  and  the 
ciorrent  drawn  through  the  two  naons  in  series  with  ths  plate  load  reaiator 
and  the  100  K.  to  -  300  v.  develops  a  neon  switching  supply  voltage  of 
~  -  100  V.  on  line  C. 

Since  line  1  is  in  the  up  state  the  cathode  of  V„  is  held  up  and 
IKGp_  is  inlriibited;  the  ca.thode  of  V^  hei.ng  up  for  the  same  reason  in- 
hibits IBM^,  and  thus  normal  drum,  activity  is  prevented  for  the  duration 
of  the  graphing  cycle.  The  choice  of  selecting  drum  heads  1  or  0  is  made 
by  the  svritch  at  the  plate  of  the  V   diode,  a  +  or  -  voltage  on  this 

'i/o   ■  '  ■  ■ 

prevented  from  being  active  during  other  than  graphing  times  because  the 
positive  controlling  voltage  for  the  l/o  line  is  derived  from  the  #1  main 
control  line,  up  only  for  graphing.  Sinrre  no  negative  voltage  on  the 
plate  of  the  V   diode  can  affect  the  l/o  line,  we  are  assured  that  the 
l/O  function  is  under  normal  operating  control  during  computing, 

A  timing  flow  diagram  ia  given  in  fig-are  ^3)  to  facilitate  follow- 
ing the  remainder  of  the  discussion  of  the  logical  process  of  the  control. 

With  T^  on  the  left  grid  of  V^n  is  negative.  Furthjer,  the  grid 

of  Y^„   (common  cathode  with  V.^  o)  is  down  due  to  control  line  1  being  down. 
Therefore  the  first  sjtic  2  pulse  will  allow  the  common  cathode  level  to  go 
negative.  This  is  used  to  ti;trn  T^^  on  (via  diode  V^-)  and  to  create  a 
false  B_  signal,  (since  the  normal  one  created  in  the  machine  is  not  now 
available)  which  is  necessary  to  turn  T„_ /-,  on  in  the  drum. 

The  false  B_  signal  is  created  in  the  following  manner.  The  left 
grid  of  V   is  down  for  the  length  of  time  of  the  sync  2  signals.  The 
right  grid  has  been  down  while  T^„  was  off,  but  now  will  rise  at  a  rate 
due  to  the  charging  of  the  EC  connected  to  it.  During  this  rise  time 
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several  "False  B^"  signals  will  hs,ve  time  to  be  initiated.  This  is  shown 
graphically  in  figure  {2k) .  Note  that  more  than  one  False  B_  signal  is 
not  harmful  if  they  occur  in  a  short  enough  interval. 

Now  that  '^■DTnnxr  ^^   ^^  ( Which  implies  that  the  counter  has  already 
done  its  first  job  of  recognizing  the  starting  block  number)  the  counter 
circuitry  is  free  to  be  cleared  in  order  to  start  its  second  Job  of  count- 
ing the  number  of  blocks  to  be  handled.  The  coincidence  circuit  input 
moat  be  looking  at  the  right  half  of  the  priming  register  in  the  drum  for 
this  second  count.  The  False  B-  signal  is  used  to  clear  the  counter  and 
also  to  turn  on  T^- /^  which  in  turn  controla  which  half  of  the  priming 
register  is  inspected. 

During  the  time  that  T^j.   is  ON  the  actual  graph  is  being  drawn 
on  the  cathode  ray  tube.  This  condition  will  continue  until  the  second 
coincidence  signal  is  obtained  indicating  that  the  proper  number  of  blocks 
of  information  has  been  handled.  On  the  receipt  of  this  second  counter 
satisfy  signal  1^^^--   is  turned  off.  The  processes  necessary  for  conclud- 
ing  the  ^aphing  cycle  and  preparing  for  the  next  cycle  are  now  initiated. 

When  Tpp  was  turned  on,  the  cathode  of  V^^-  sent  a  down-going  sig- 
nal* to  the  drum's  null  order  gate  output  line,  inhibiting  any  null  order 
gate  signal  from  being  transmitted  to  the  main  machine.  With  T„^-_j-  OFF 
(at  second  counter  satisfy  time)  in  coincidence  with  T„_ /^  being  ON  a 
positive  null  order  gate  signal  is  obtained  froa  the  drum  circuitry. 
This  causes  Vq  to  conduct  and  turn  cf f  T^^ .  The  T„-_  control  tube  in 
the  drum  is  then  released  (via  V^,)  so  that  Tg^  is  turned  OFF.  With  Tg_ 
OFF,  T^- /_  is  turned  off  and  as  a  consequence  the  null  order  gate  request 
signal  is  turned  off.  At  this  time  T__  is  turned  off  via  tubes  V_  and  V-_. 
Note  that  it  was  necessary  that  T-_  remain  on  until  the  previous  null 
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order  gate  signal  sequence  had  finished.  Tliat  Is,  the  NOG  inhibition  is 
used  to  cover  the  entire  turning  off  period.  This  completes  the  graphing 
cycle  and  all  elements  of  the  control  have  "been  restored  to  normal  or 
standby  condition. 

Before  starting  the  next  cycle  certain  conditions  must  be  met. 
These  are  incorporated  in  tubes  V,,  and  V  .  If  an  IWG  signal  has  been 
waiting,  then  T^,  is  prevented  from  coming  back  on  to  start  a  new  graph- 
ing cycle.  In  this  case,  INGp  will  initiate  the  appropriate  external 
order.  Note  that  if  an  ING  request  comes  during  a  graphing  cycle,  the 
computer  will  wait  until  that  cycle  is  completed  before  the  order  is 
executed. 

A  new  graphing  cycle  may  not  begin  until  the  completion  of  the 

external  order  (i.e.  end  of  the  null  order  gate  signal).  At  this  time 

the  common  cathodes  of  V^ ^  and  V^  will  become  negative  turning  T^,  on 

J--L       lii  Gl 

via  Vpj^,  initiating  a  new  graphing  cycle.  Note  then  that  graphing  op- 
erates continually  except  for  interruptions  by  external  machine  orders, 
that  is  drum  or  IBM  operations. 

The  remaining  functions  of  the  graphing  control  are  to  provide 
proper  signals  for  the  deflection  and  the  beam  turnon  routine.  The  line 
marked  "Graph  Routine",  connected  to  the  cathodes  of  V^  and  V^o,  i3 
used  for  this  purpose.   It  was  seen  that  this  line  will  be  down  for  the 
simultaneous  conditions  of: 

1.  Graphing  cycle  in  progress  (TG,  on) , 

2.  T___^^  ON  (information  being  read  out  of  the  drum) ,  and 

BLOCK 

3.  Sjac  2   pulses  being  transmitted  (timing  pulses  for  each 
word  of  information) . 

Thus  a  beam  turnon  signal  is  available  for  turning  on  the  cathode  ray 
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tube  during  each  desired  word  readout. 

The  Graph  Eoutine  line  is  also  used  to  provide  the  clear  and  gate 
signal  for  the  deflection  input  circuitry.  This  signal  (during  graphing 
replacing  the  computer  clear  signal  used  during  slave  operation)  acts 
as  a  trigger  for  the  clear  driver  of  the  gi-aphing  deflection  circuit. 
Since  each  clear  and  gate  action  inserts  nev  digit  information  into 
the  deflection  adder,  then  the  total  result  is  that  at  the  time  of  each 
vord  read  out  for  graphing  the  transfer  of  that  vord  into  the  graphing 
deflection  circuit  is  accomplished.  Beam  turnon  at  that  time  then  dis- 
plasrs  the  information  in  that  word  as  a  point  at  the  desired  location  on 
the  screen  of  the  cathode  ray  tube.  Subsequent  words  then  extend  this 
point  to  a  line  or  lines,  and  graphing  is  accomplished, 
e.  Division  spill  check. 

Proper  machine  numbers  m-ast  lie  in  the  range 
-1  <  X  <  1  and  the  division  circuitry  will  give  a  grossly  wrong  answer 
if  a  proposed  quotient  q  =  x/y  lies  outside  this  range.  A  check  circuit 
was  added  to  stop  the  machine  if  e:a   improper  division  is  attempted. 

The  division  process  used  consists  of  39  similar  steps  determin- 
ing a  sign  and  38  quotient  digits.  It  is  a  restoring  method  in  which 
the  adder  output  is  examined  for  overdraft  before  it  replaces  the  old 
remainder.  In  case  of  overdraft,  the  adder  output  is  not  used  and  the 
old  remainder  is  doubled.  The  first  step  of  a  proper  division  should 
always  present  an  overdraft,  and  this  fact  is  used  for  the  spill  check. 
(See  figure  (26)) . 

As  the  division  order  is  entered,  the  initial  memory  sequence 

which  brings  the  denominator  into  R  sets  the  division  spill  toggle 

2 
to  "1"  via  Eed  Clear  R  .  Note  that  this  happens  for  all  memory  orders 


Figure  26  -  1-1*9. 
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but  that  the  right  hand  gate  tube  will  immediately  unset  the  toggle 
unless  a  division  order  is  present.  The  "1"  state  of  this  toggle  is 
used  to  inhibit  the  "accept"  (i.e.  adder  path)  clear  and  will  stop  the 
machine  if  such  a  clear  is  requested.   In  a  proper  division,  however, 
the  first  step  is  a  "reject"  (due  to  the  overdraft)  and  the  "reject" 
clear  -Black  Clear  E  -  resets  the  toggle. 

At  the  time  of  a  division  spill  stop,  the  original  numerator  is 

3 
still  in  E^  and  the  denominator  in  K  so  that  the  operator  can  judge 

which  is  out  of  scale.  The  spill  inhibition  is  removed  by  turning  off 

the  order  gates  (the  standby  state  of  the  control  box) . 

f .  Zero  shift  order . 

In  the  initial  design  of  the  machine  it  was  assumed 
that  if  a  shift  order  was  entered  at  least  one  shift  would  be  done  and 
so  the  check  for  shift  count  satisfaction  was  placed  at  the  end  of  the 
shift  cycle.  For  a  number  of  reasons  a  zero  shift  was  requested  by  the 
users  but  primarily  for  routines  in  which  the  shift  argument  is  varied 
by  the  code. 

The  zero  shift  was  accomplished  by  moving  the  sensing  for  count 
satisfaction  to  the  beginning  of  the  shift  cycle  for  all  except  input- 
output  orders.  When  executed,  a  zero  shift  order  leaves  the  contents 
of  both  E^  and  E  undisturbed. 

g.  Memory  high  voltage  change. 

The  IAS  memory  was  designed  for  an  accelerating 
voltage  of  1000  volts  on  both  second  and  third  anodes,  this  being 
originally  a  compromise  between  spot  definition  (which  improves  with 
increasing  voltage)  and  storage  flaws  (which  become  worse) .  As  experi- 
ence was  accumulated,  both  in  the  machine  and  in  the  CRT  test  rack,  it 
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became  evident  that  the  compromise  voltage  should  have  been  set  higher: 
flaws  have  not  been  a  problem  but  read-around  has.  This  is  true  partly 
because  flaws  can  be  completely  evaluated  in  the  test  rack  whereas  read- 
around  cannot;  and  partly  because  a  flaw-free  tube  stays  so  throughout 
its  life  while  read-around  slowly  deteriorates. 

Unfortunately  both  grid  and  deflection  plate  circuits  are  direct 
coupled  in  the  memory  and  any  significant  change  in  cathode  to  second 
anode  voltage  involved  a  circuit  change  rather  than  merely  a  power  supply 
change.  However,  preparations  were  made  during  the  machine  move  (see 
C,l.,  a.,   i.)  for  a  higher  voltage,  and  all  incoming  cathode  ray  tubes 

were  flaw  tested  at  2000  volts  E,  rather  than  1000  volts. 

b 

In  addition  since  the  old  high  voltage  supply  was  not  capable  of 
delivering  a  higher  voltage,  a  new  supply  was  designed  and  built  (see 
figure  (27) ) •   It  does  not  differ  in  any  important  respect  from  the  old 
one  except  for  the  higher  voltage,  a  more  easily  serviced  layout,  and  an 
ungrounded  negative  retiirn.  Of  some  interest  may  be  the  use  of  the  gas 
filled  3B2U  rectifiers.  These  have  proved  much  more  reliable  than  the 
8016 's  previously  used  and  no  trouble  is  caused  by  the  ionization  transi- 
ents of  the  gas . 


Figure  27  -  l-kS-. 
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REGULATED  VOLTAGES 


Nominal 

Adjustment 
Range 

Set  At 

M/C  Min.   ma. 

m/C  Max.    ma. 

Normal  No- 
Load  Volt 

Ref.   Volts 

+  86 

80-  95 

86.0 

-   170 

10 

+100 

+  82 

+100 

96-108 

100.0 

0 

10 

+138 

+  82 

+120 

111-131 

125.0 

0 

10 

+155 

+  82 

+136 

12ii-li^8 

136.2 

-   200 

50 

+152 

+  82 

+150  D 

l!;0-l6l 

151.0 

810 

1160 

+2^0 

+  82 

+150  G 

■Lk3-l6^ 

150.0 

150 

510 

+205 

+  82 

+150  M 

ii^o-i6o 

152.0 

870 

1120 

+2^6 

+  82 

+161 

153-168 

163.0 

-      35 

275 

+230 

+  82 

+186 

1714.-205 

19^^.1 

-        6 

20 

+239 

+  82 

+211 

203-216 

210.8 

0 

150 

+250 

+161^ 

+220 

215-228 

223.0 

510 

610 

+3iv5 

+l6i^ 

+300  A 

293-313 

303.5 

iv5 

J^35 

+376 

+161^ 

+300  D 

292-316 

303.0 

37 

1^5 

+362 

+16!+ 

+300  V 

290-312 

293-5 

lll-OO 

1980 

+369 

+l6i+ 

+336 

330.31^1+ 

33^^.0 

700 

1150 

+380 

+16^ 

-17^ 

167-177 

170.0 

595 

650 

-238 

(-300) 

Nominal 

Load  and 

No-Load 

Volts 

REFERENCE 
"No-Load"  MA 

]  AND  SCREEN  VOLT 
"Loaded"  MA 

AGES 

Deaigned 
Min.  Ma. 

Deaigned 
Max.  Ma. 

+  82 
+l6i^ 
+232 
+3lif 

81.9 
I6I1.2 
229-0 
312.0 

2 
3 
0 
0 

19.0 

10.0 

5-5 

3.5 

10 
6 
3 
2 

28 

15 

11 

6 

Regulated  Supplfca 

27  -  6080 

2  -  6336   (ir 
16  -  6au6 

(16  in  all) 
I  "300  v") 

Total  In 

Ref.  &  Screen  Volts 

Cliassis 

r 
C 

r 
C 

c 
r 
C 

ternal  Watts 

I  -  6au6 

;   -  12BH7 
I  -  0A2 
I  -  5651 

Filam 
D.CD 

ents     560 
iss.   1100 

Figuxe  28. 
Table   of  Regulatad  Power  Supplies. 
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3-  Machine  operation. 

a.  During  the  year  covered  by  this  report,  the  machine 
was  operated  on  a  fully  scheduled  basis  except  for  the  four  month  period 
in  which  it  was  moved.  During  the  eight  active  months,  the  machine  was 
run  for  a  total  of  326^4-  hours;  of  which  10^  was  spent  in  routine  mainte- 
nance, 17%   in  engineering  improvements,  9^  in  unscheduled  maintenance, 
and  6ki^   in  productive  operation. 

The  overall  machine  efficiency  for  the  year  could  be  expressed  as 
the  productive  operation  time  divided  by  total  "on"  time  with  engineering 
excluded.  As  a  percentage,  this  "economic"  efficiency  was  77^-   In  judg- 
ing the  state  of  engineering  control  of  the  machine,  it  was  perhaps  more 
meaningful  to  exclude  also  the  routine  maintenance  time  as  well  as  engi- 
neering improvement  time  from  the  denominator.   If  this  is  done,  we  obtain 
an  "engineering"  efficiency  of  87^  for  the  year. 

b.  Weekly  time  records  for  the  year  are  given  in  Table  I. 
For  completeness  the  time  records  for  the  previous  year  are  given  in 
Table  II.  Eesults  for  the  two  years  are  summarized  in  Table  III,  by 
quarters,  the  last  quarter  of  1952  being  the  first  for  which  these 
records  were  kept.  With  respect  to  all  tables,  the  time  breakdowns  are 
defined  as  follows: 

Eoutine:  All  time  spent  performing  a  standard  set  of  tests  pre- 
liminary to  each  day's  operation. 

Engineering;  All  time  spent  in  engineering  improvements  or 
additions  to  the  machine . 

Unscheduled  Maintenance ;  All  time  spent  in  diagnosing  and  re- 
pairing conditions  which  have  interrupted  computation  or  which  threaten 
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to  cause  errors.  Alao  included  ia  bad  operating  time  which  is  clearly- 
due  to  machine  trouble. 

Operation;  All  machine  "on"  time  not  accounted  for  by  the  above. 


Table  I. 

I- 

-56. 

BASIC  EECORDS 

DERIVED 

Week  Of 

Routine 

Engineering    Unscheduled 
Maintenance 

Operation 

Total 

(UnM  +  Op) 
Available 

(Op  /  Av) 
^  Operation 

29  June 

7.2 

9.3                     11.3 

125.2 

153.0 

136.5 

92 

6  July 

12.9 

3.5                     19.8 

99.8 

136.0 

119.6 

83 

13  July 

13.3 

8.6                       8.5 

103.8 

13*^.2 

112.3 

92 

20  July 

8.7 

3.h                    11.7 

I2I+.2 

11^8.0 

135.9 

91 

27  July 

8.0 

0                    3.0 

57.2 

68.2 

60.2 

95 

3rd  Quarter 

50.1 

2lf.8                   5ii..3 

510.2 

639'k 

561^.5 

91 

Machine  off  for  moving  1  August  1953  to 

10  November  1953. 

10  Nov. 

0 

32.0                   0 

0 

32.0 

0 

_  « 

l6  Nov. 

0 

38.7                  0 

17.1 

55.8 

17.1 

_  * 

23  Nov. 

0 

6k. 0                   0 

0 

6k.O 

0 

_  « 

30  Nov. 

0 

67.1                   0 

12.9 

80.0 

12.9 

_  * 

7  Dec. 

114-.9 

13.2                    7.0 

k3.6 

78.7 

50.6 

86 

ll;  Dec. 

11.9 

15.5                    6.8 

1^5-8 

80.0 

52.6 

87 

21  Dec. 

8.i^ 

2.0                    2.7 

27.8 

i^O.9 

30.5 

91 

It-th  Quarter 

35-2 

232.5                   16.5 

11^7-2 

U31.i+ 

163.7 

90 

k  Jan. 

9.7 

9.1                  23.3 

1^0.5 

82.6 

63.8 

63 

11  Jan. 

lit.3 

1.0                   13-0 

1^7.5 

75.8 

60.5 

78 

18  Jan. 

10.1 

2.6                       li^.2 

63.1 

90.0 

77.3 

81 

2U  Jan. 

8.9 

16.6                   0 

65.6 

91.1 

65.6 

100 

1  Feb. 

7.0 

31.7                   ^'9 

U0.8 

eh.k 

lv5.7 

89 

8  Feb. 

9.7 

8.7                    2.3 

63.9 

eh,6 

66.2 

96 

15. Feb. 

ll.l^ 

19.1                    3.7 

I13.8 

78.0 

1^7.5 

92 

22  Feb. 

6.2 

8.3                    2.6 

59.9 

77.0 

62.5 

96 

1  Mar. 

13.3 

13.2                    2.0 

52.6 

81.1 

5U.6 

96 

8  Mar. 

11.2 

12.7                  11-1 

51^.7 

89.7 

65.8 

83 

15  Mar. 

7-3 

8.0                    5.3 

68.2 

88.8 

73.5 

93 

22  Mar. 

6.5 
115.6 

7.3                    3.3 

65.8 

82.9 

69.1 

95 

1st  Quarter 

138.3                  85.7 

666. k 

1006.0 

752.1 

88 

29  Mar. 

6.k 

0                     27.3 

69.2 

102.9 

96.5 

72 

5  Apr. 

10.1 

5.7                    1'^ 

57.7 

80.9 

65.1 

88 

12  Apr. 

11.5 

808                    k.7 

U6.1 

71.1 

50.8 

91 

19  Apr. 

10. 1^ 

17.7                       20.0 

38.3 

86.4 

58.3 

66 

26  Apr. 

10.5 

11.7                       22.3 

1^7.2 

91-7 

69.5 

68 

3  May 

8.8 

11.3                         0.7 

69.5 

90.3 

70.2 

99 

10  May 

9.7 

27.0                       12.3 

50.3 

99.3 

62.6 

80 

17  May. 

8.9 

17.0                         3.1 

62.8 

91.8 

65.9 

95 

2k  May 

9.9 

16.2                         7.1 

55.0 

88.2 

62.1 

88 

31  May 

11.1 

0.5                        10.2 

55.2 

77.0 

63.k 

8k 

7  June 

5.5 

32.8                        16.2 

57A 

111.9 

73.6 

78 

ik  June 

Ik.l 

1^.5                       10.1 

73.2 

102.5 

83.3 

88 

21  June 

7-1 

3.2                         9.7 

72.8 

92.8 

82.5 

83 

2nd  Quarter 

121^.6 

156. 1^                     151.1 

75i^.7 

1186.8 

905.8 

83 

djustment  period  after  move. 
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Table 

1  n. 

1 

BASIC  RECOBDS 

DERIVED 

Week  Ending 

Eoutine 

Engineering 

Unscheduled 

Operation 

Total 

(UnM  +  Op) 

(Op  /  Av) 

Maintenance 

Available 

^  Operation 

31  Oct. 

14.5 

32.5 

6.0 

26.0 

79-0 

32.0 

81 

7  Nov. 

11.0 

22.5 

6.0 

V1..0 

80.5 

1+7.0 

87 

Ik  Nov. 

10.5 

17.0 

37.5 

liv.5 

79.5 

52.0 

28 

21  Nov. 

9.5 

2V.0 

27.5 

18.0 

79-0 

1^5.5 

1+0 

28  Nov. 

9-5 

11.5 

10.0 

33.0 

61+. 0 

1^3.0 

77 

5  Dec. 

10.5 

8.5 

9.5 

51.5 

80.0 

61.0 

81+ 

12  Dec. 

10.0 

19.0 

9.0 

J4-2.0 

80.0 

51.0 

82 

19  Dec. 

11.0 

27-5 

5.5 

36.0 

80.0 

1+1.5 

87 

26  Dec. 

6.5 

5.5 

1.5 

23.0 

36.5 

2I+.5 

9k 

2  Jan. 

- 

-     (Holiday)    ~ 

- 

- 

- 

Total 

93.0 

168.0 

112.5 

2B5.O 

658.5 

397.5 

73 

9  Jan. 

9.2 

18.1 

26.9 

25.8 

80.0 

52.7 

k9 

16  Jan. 

8.3 

1^.8 

35.7 

31.2 

8c. 0 

66.9 

hi 

23  Jan. 

10.1 

22.5 

21.8 

25.6 

80.0 

kj.k 

3h 

30  Jan. 

9.5 

37-0 

17.6 

17.7 

81.3 

35.3 

50 

6  Feb. 

8.5 

36.3 

6.8 

28.5 

80.1 

35.3 

81 

13  Feb. 

7.8 

29.7 

9.3 

3*^.9 

81.7 

1+1+.2 

79 

20  Feb. 

7.8 

3.9 

6.1 

62.2 

80.0 

68.3 

91 

27  Feb. 

8.l^ 

2.5 

10.7 

58.1 

79.7 

68.8 

81+ 

6  Mar. 

8.1 

0 

18.6 

65.6 

92.3 

81+. 2 

78 

13  Mar. 

10.5 

26.0 

10.2 

57.3 

lOl+.O 

67.5 

85 

20  Mar. 

12.5 

6.1 

0.1 

85.1 

103.8 

85.2 

99 

27  Mar. 

7.8 

10.2 

9.2 

77.0 

101+.2 

86.2 

89 

Total 

108.5 

197.1 

173.0 

569.0 

lOi+7.6 

7'+2.0 

71+ 

2  Apr. 

I'k 

6.k 

k.6 

69. J^- 

87.8 

74.0 

91^ 

12  Apr. 

13.1 

16.1 

1.  n 
-<■•  1 

78.7 

112.6 

83.1+ 

9h 

17  Apr. 

9.1 

8.2 

21.5 

65.1 

103.9 

86.6 

75 

2k  Apr. 

10.8 

7.7 

3.3 

76.0 

97.8 

79.3 

96 

1  May 

9.2 

11.8 

3.3 

37.7 

112.0 

91.0 

95 

8  May 

7.3 

11. U 

39.1 

1+4.1 

101.9 

83.2 

53 

15  May 

9.2 

2.k 

kk.^ 

k'J.9 

104.0 

92.1+ 

52 

22  May 

11.9 

0 

11.6 

80.5 

lOl+.O 

92.1 

87 

29  May 

12.7 

2.0 

32.9 

83.6 

131.2 

116.5 

72 

7  June 

8.3 

8.!+ 

7.8 

82.6 

107.1 

90.1+ 

91 

Ik  June 

12.3 

15.7 

23.9 

97-3 

1^9.2 

121.2 

80 

21  June 

11.8 

9'k 

5.iv 

105.9 

132.5 

111.3 

95 

28  June 

10.8 

1-h 

5.5 

12^.0 

1^7.7 

129.5 

96 

133.9 


106.9 


208.1 


101+2.8 


11+91.7 


1250.9 


83 
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Table   III. 

BASIC 

EECOEDS 

DEEIVED 

Quarter 

Eoutine 

Eng. 

Unsched. 

Operation 

Total 

(UnM  40P) 

(Op 

/Av) 

Maint . 

Available 

A 

Op 

eration 

1952 

i^th 

93.0 

163.0 

112.5 

285.0 

658.5 

397.5 

73 

1953 

1st 

108.5 

197-1 

173.0 

569.0 

10i^7.6 

7^+2.0 

Ih 

2nd 

133.9 

106.9 

208.1 

lOi+2 . 8 

lli91.7 

1250 . 9 

83 

3rd 

50.1 

2i|.8 

^h.3 

510.2 

639-^ 

561^.5 

91 

i^th 

35-2 

232.5 

16.5 

II+7.2 

i^31.i^ 

163-7 

90 

195^^ 

1st 

115.6 

138.3 

85.7 

666. k 

1006.0 

752.1 

88 

2nd 

12i^.6 

156.!+ 

151.1 

T?h.7 

1186.8 

905.  a 

83 

1953-5^  totals    (3rd  and  Uth  quarters   of  1953;    Ist  and  2nd   quarters   of  19^h) 
325.5  552.0         307-6         2078.5  3263.6       2386.1 


10^ 
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PART  II  -  MATHEMATICS 


II-l. 


A.     BLAST  WAVE  CALCULATION* 

0.  Introduction. 

In  the  succeeding  pages  we  discuss  both  the  numerical 
analysis  and  the  calculations  we  performed  in  connection  with  a  blast 
wave  computation.  The  problem  itself  is  that  of  studying  the  decay  and 
propagation  of  the  shock  wave  generated  by  a  very  strong  point-source 
explosion  in  an  ideal  gas  under  the  assumption  of  spherical  symnetry. 
We  additionally  were  interested  in  the  behavior  of  the  pressure  as  a 
functional  of  radial  distance  for  fixed  times  and  of  pressure  as  a 
function  of  time  for  fixed  distances . 

The  method  employed  for  handling  the  conditions  of  Eankine  and 
Hugoniot  is  an  iterative  one  and  is  discussed  in  section  3  below.  Re- 
cently Erode  of  the  Rand  Corporation  undertook  a  blast  wave  calcula- 
tion  closely  paralleling  the  one  described  in  this  paper."  His  method 
is,  however,  quite  different. 

In  the  succeeding  sections  we  discuss  the  physical  problem,  the 
numerical  algorithm  and  the  computational  results.  The  latter  are  con- 
tained principally  in  a  set  of  graphs  appended  at  the  end  of  the  paper. 

1.  The  differential  equations. 

It  is  convenient  to  let  Po     ,     po       ^e  "ti^e  original  den- 
sity and  pressure,  i.e.  the  density  and  pressure  in  the  undisturbed  gas; 
t  is  time,  r  is  the  Lagrangean  radius,  R  the  Eulerian  radius,  P   the 
density,  p  the  pressure,  Jl   the  position  of  the  shock,  D  the  velocity 


*  This  paper  has  been  submitted  to  a  mathematical  journal  for  publication. 

#  Erode,  H.  L. ,  Results  of  a  Numerical  Integration  of  a  Spherical  Shock 
from  a  Point  Source,  Rand  Report.  (This  paper  contains  references  to 
the  literature  on  the  subject. 
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of  the  ahock,  and  U  the  mass-velocity  immediately  back  of  the  shock- 
Then  Newton's  second  law  implies  at  once  that  behind  the  shock 


The  equation  of  continuity,  i.e.  of  the  conservation  of  mass,  requires 
that  everywhere 

Finally  the  equation  of  state  requires  that 

P  =  ^  f         > 
where  k  is  a  function  of  r  but  not  of  t  behind  the  shock,  i.e. 

k  =  k(r). 
These  equations  describe  the  conditions  behind  the  shock  but  do 
not,  as  we  know,  hold  across  the  shock.  Instead  we  must  consider  the 
well-known  equations  of  Eankine  and  Hugoniot.  They  are  these: 


u  h4    '^-^ 


-^°    -/ff*ij^-^/'r-/)  po 


P     =.A 


^^^ff^j)p^fr-op." 


The  boundary  conditions  are  two  in  number,   one  being  at  the  origin, 
i.e.   at  r  =  0  and  the  other  at  the  shock,    i.e.   at  r  =    it    •     They  are 

E(0,  t)    =  0  for  all  t, 

E(    ^    ,  t)    =  ^  for  all  t. 
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where  of  course,  Jt  is  a  function  of  t, 

The  initial  conditions  are  discussed  in  §!<•  below. 

In  passing  we  note  that  the  equations  (1.5),  (1.6),  (1.7),  (1-9) 
are  not  independent.  If  we  differentiate  (1.9)  with  respect  to  t,  we 
find 


2 


di?     ^     9e       d^ 


dr         dt  dt         dt 

But  hj  (1.2)   and   (19) 

Pe 


/■f!  =  /c 


at  the  shock.     Thxia    (1.11)    is  equivalent  to 

(1  -  ^)   D  =  U. 

It  ia  trivial  now  to  see  that  this  relation  plus  any  two  of  (1.5),  (1.6), 
(1.7)  implies  the  third. 

2.  The  difference  eqixations. 

In  what  follows  we  shall  show  how  we  replace  our  original 
implicitly  formulated  problem,  equations  (1.1)  -  (1.9)^  ^7  a  finitiatic, 
explicit  formulation  which  approximates  the  original  one. 

We  consider  first  (1.1),  (1.2).  We  combine  them  into 

P)i^     ~  p  dn     dr  yj  r"-    dr 

Thus  our  sjrstem  (1.1)  -  (I.9)  becomes  now  (2.1),  (1.2)  -  (1.9). 

As  will  be  seen  later  it  is  convenient  to  introduce  a  change  of 
▼ariables.  This  change  is  the  following: 

8  =  r2/2 
In  terms  of  this  new  variable  (2.1)  becomes 
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d'R     ^    _    J_       7^^  dp 

■where  we  now  understand  E,  p  to  be  functions  of  s,  t.  We  also  use  the 
quantity 

in  what  follows. 

We  now  replace  our  continuous  variables  r,  s,  t  by  discrete  ones. 
The  time  variable  t  is  replaced  by 

t^    (h  =  0,  1,  ...) 
where 

t^^^  =  t^+  At''^^'^, 

with  the  A  t^"^  ■'■'^  determined  as  indicated  in  (3  =  35)  below.  We  note 
that  the  t^  (h  =  0,  1,  ...)  are  not  necessarily  equidistant.  The  spatial 
variables  are  defined  for  each  t^  in  this  fashion.  The  distance  from  0 
to  y^^  =  j/^(t^)  is  divided  into  L  equal  parts,  thus; 

ds^  =  y  /l     (L  a  positive  integer) . 
Then  a .  is  given  with  the  help  of  the  relation 
3  .  =  jds      (j  =0,  1,  ... ,  L) ; 
note  that  the  s .  are  functions  both  of  j  and  h.  We  suppress  the  depend- 

J 

ence  on  h  when  no  ambiguity  can  arise. 

By  this  means  our  spatial  scales  cover  just  the  relevant  portions 
of  space,  namely,  a  gradually  expanding  region  about  the  center.  In 
order  to  compare  results  at  two  different  time  levels  it  is,  of  course, 
necessary  to  take  into  account  the  dependence  of  the  s^  on  h. 

These  things  being  understood  we  discuss  the  difference  equations 
to  be  used  in  place  of  (2.3),  (1-2)  -  (1.9)- 
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Consider  the  formula 

?     '^1    *j2ij(^l-i-P}*'4)^     ds'^  2  els' 

(j  =1,  2,  .o.,  L-1). 

In  this  R.    corresponds  to  E.    except  that  the  readjustment  of  the 

spatial  scale  to  the  change  from  v   to  :7     is  not  contained  in  it, 

i.e.  it  corresponds  to  the  same  Lagrangean  positiono  To  understand  this 

recall  that 

:^  h+1   ^/ . ,  h  .h+l» 
Bj    =  E(jds  ,  t   ) 

whereas 

_  h+1   _,  . ,  h+1  .h+lv 
E.    =  E(jds   ,  t   ). 

Further  E  .^"'"■''  is  the  estimate  of  R .  "*"  obtained  by  extrapolating  linearly 
J  J 

in  t  from  E.   ,  E.  ,  taking  all  three  quantities  in  the  spatial  scale 
at  time  t  .   If 

~~  h-1   ^/ . ,  h  .h-l» 
E.    =  E(jds  ,  t   ), 


then 


Kj      "  ^^  -^  ^^  k-'/z    (^1  -l<i    J  ,  (j  =1,  2,  ...,  L-l) 


It  is  thus  clear  that 


di 
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We  consider  nert  the  remaining  terms  in  (2.9)-  With  the  help  of  (2A) 
and  (2.7)  we  find  that 


Jzl/^n-^^     Pf^^l    ^/.      ~~d^  2ds^ 

Thus  through  second  order  terms  formula  (2-9)  is  an  approximation  to  (2.3) • 

We  turn  attention  now  to  relations  (1.2),  (l«3)-  Since  k  is  in- 
dependent of  t  ,  we  see  that 


where  the  bars  over  p  and  f    Imve  the  significance  explained  earlier  in 
connection  with  E,  i.e. 

We  find  with  the  help  of  (1=2)  that  (2.11)  becomes 


(j  =  0,  1,  ...,  L-l) 


The  remaining  relations  (1-5)  -  (1-9)  are  already  finitistic  in  character 
and  require  no  emendation. 

Our  new  system  is  thus  (2.9),  (2-15),  (1-5)  -  (1»9)  together  with 
the  definitions  (2.10),  (2.11),  (2.12),  (2.li^).   It  remains  only  to 
describe  the  computational  algorithm  based  on  these  equs-tions  and  the 
process  used  for  obtaining  the  initial  values. 

3.  The  computational  algorithm. 

We  assume  as  given  the  following  quantities; 
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Ej^,   Pj-l/2   iJ   -  0,   1,    ...,   L),         B^.^""^   (j   =  1,   2,    ...,   L-1); 

^^        D  ^-'/^     D  ^-3/2,   t\   dt^^^/2,   dt^-3/2   „ 

We  describe  in  this   section  how  to  calculate  the  corresponding  quantities 

.    ,  .        . h+1 
at  time  t 

First  we  notice  that  formula   {2-9)   enables  us   immediately  to 

define 

E.^-'^  (j  =  1,   2,    ...,  L-1); 

(1.8)    yields 

■  bV^^  .  0. 

Thus   only  E  is   lacking  --  this  value  will  be  determined  with  the 

Li 

help  of  the,  as  yet  unused,   relations    (1.5)    (l-7)>    (l°9)<' 
Next  formula   (2.15)    permits  us  to  calculate  easily 

PjW2  (^    =°'    "■'    —^-2) 

from  the  sequence    (3.3) >    (3°^)   above^      Given  E  we   can  also  form 

-  h+1 
Pl-1/2 

directly  from  (2.15).  Finally  by  the  spherical  symmetry  of  the  problem 

-  h+1   -  h+1 
P-1/2  =  P  1/2   ° 

We  turn  attention  now  to  the  determination  of  E^^   ,  i.e.  to  the 
utilization  of  the  shock  conditions  (1.5)  -  (1-7),  (1.9)°  The  procedure 
to  be  described  is  an  inductive  one,  as  will  be  seen.  It  is  due  to  E. 
Peierls  originally.  The  induction  insofar  as  the  shock  velocity  is  con- 
cenied  yields  a  value  ^  -,'^   for  this  velocity  given  ^o^    and  y^_-^     • 
Thus  initially  one  requires  estimates  .  o^  ,  ,  o<y  « 
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It  starts  in  this  manner:  Given  the  quantities  (3'1)>  (3-2)  it 
is  desired  to  find  first  approximations  oO    ,   ,  oU    to  D   '  •  This  is 
done  with  the  help  of  the  relations 

We  thus  have  first  approximations   oZ?  ,  o^y  to  the  shock  velocity  at 

^h+l/2   ,h  .  ,.h+l/2 
time  t   '   =  t  +  dt   '  . 


The  general  step  in  the  induction  can  now  be  described.  Given 
an  approximation   o3  to  D  ^     we  find  from  (1-7)  a  value  p  of  the 
shock  pressure  at  t   '  j  it  is  this 

P  =   '     Po  D^-  ^-^   Pp. 

Given  this  estimate  for  the  shock  pressure  (1.6)  combined  with  (l.?) 
yields  a  value  ,U  of  the  mass-velocity  at  the  shock;  it  is  this 


P  -  P, 


U  = 


0 


/o 


D 


Finally  relation  (1-12)    implies  that 

fo     ,   u 

Thus  an  estimate   P   for  the  shock  density  can  be  obtained. 
We  describe  now  how  these  quantities 

k^^  kP'  k^^  kf  ' 

are  used.  From  what  has  preceded  we  know  the  values  of 

^  h+1  =•  h+1 
^L-2  '  \-l 

and  thus  we  know  , 
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i.e.  we  have  U,  the  mass  velocitj,  at  (L-2)  ds  ,  (L-1)  ds'^  and  at  the 
shock  u     ^     ' ".  From  these  we  ;an  by  quadratic  interpolation  find  U  at 
Lds  ■ .  This  is  the  piorpose  of  the  formula 

(Notice  in  theae  formulas  that  T    ,     -U    ,     a    and  thus  Sr   depend  upon 

li 

the  inductive  index  k.)  With  the  help  of  .  O  t  we  now  define  an  estimate 

K   L 

for  B_    by  means  of  the  formula 

Li 

^  h+1       ^    h  r 

Wi1;h  this  quantity  ws  can  form 

-  h+1 
k^L-l/2 

from  (2.15) •  But  we  already  know 

-  h+1     -  h+1 

thua  by  a  quadratic  extrapolation  we  can  find  an  estimate 

kP(h+l/2) 
for  the  pressure  at 

B  =  /^-^^/^  t  =  t^^l   . 

This  quantity  is  obtained  from  the  formula 
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with 


To  see  this  value  for     u     is   correct,  we  note  that  the  s -distance  in  ques- 
tion is  from  s  =   (L  -   3/2)   ds "  to  a   =  7    "^       •     This   ia 

ds^^  .  dsV2  +    f^^'--  -  L  d  s^  = 

=     ^(3   +  L   (2  y     +   ^2/2))    ds'^      . 
The  reader  is  referred  to  Figuie  I  'below. 


-» )( ' St • fr- 


L-5/2  L-3/2  L-1/2 


L"3  L-2 


L-1 


h+1 


D''"'/^  p,  ^0'  ° 


Fig'ore  I. 


we  go  now  to  obtain  another  estimation 

k^(h+l/2) 
for  the  pressure  at  the  (a,  t) -point  (3-l8).  We  do  this  with  the  help 
of  the  adiabatic  law.  The  value  of  p  in  question  may  be  obtained  as 


k^(h+l/2)   k^ 


^ 


''*')^-<.^r'' 


where 
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C  '  /o//  -  1  ■  J"'yP-  91"'''  =5?"  -  .»  •  ^*''*''"    ■ 

(Note  that  both    S     and  ^         are  functions   of  k,   the  inductive  index.) 

Given  these  two  not  necessarily  equal  values   of  p(jr  )  "^       ) 

we  are  able  to  form  an  "error  indicator" 

k^     ^  k^Cli+1/2]      '  k^(h-l/2) 
and  to  uae   it  aa  a  means   of  obtaining  a  new,    improved  value  ^.-.D  for  the 
shock  velocity.     This   is  done  with  the  help  of  the  formulae 

and  (3.8). 

The  procedure  is  iterated  until  for  some  pre-assigned  £"  we  have 

When  this  point  is  I'eached,  by  definition 

Dh+1/2__^D.  ^^^^=  ^\d^-^""/2  «  dt^-^^/2,  B^^^^=iE^^^\  t^^^=t^Vdt^-^l/2   . 

In  (iS.S),  (3.^)  "we  saw  how  B.^'^^   (j  =0,  1,  ...,  L-1)  were  obtained. 
We  now  have  the  complete  sequence 

R^.^'""^      (J  =  0,  1,  ...,  L) 


and  also  from  (3o)j  (3«6),  O-T) 


-  h+1 


P 


-1/s 


(j  =  0,  1,  ...,  L) 


We  describe  next  how  to  obtain  R.   ",  P^.^^  (J  ^  ^'  ■'■»  "  ° ' '  ^^  °^*  °'^ 
these  sequences.  This  is  done  with  the  aid  of  interpolation  formula 

(3.19)- 

To  form 

E.^^^     (j  =  1,  2,  ...,  L-l) 
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ve  use  (3 •19)  with 

f .  =  E  .^-\  u  =  j  (2  >^  ^  ^2)  ^   __  ^^^1/2   .   ,,h^l/2/^  h 

Thug  we  have  obtained,  as  desired,  the  sequence 

E.^'""'^     (j  =  0,  i,  ,..,  L). 

h+1 
Similarly  we  find  the  p.    ,  /^   .     Specifically  by  setting 

f.;  =  Pj_i/2  .   ^=  (J   -  1/2)  (2   V+   y^) 
ve  obtain  p.^^w^   (j   =  1,   2,    ...,   L-1)    from  (3-5),    (3«6),    (3.?).     Then 


h+1  h+1 

P-1/2       =     ^/2        • 

The  final  value  Pt'T/o  ^^   obtained  by  setting  u  =  1  +   (L  -  l/2)(2-y'+    y>   ), 

i.e.   by  an  ertrapolatiorio      Thus  ve   obtain 

PJ-V2  0=0,    1,    ...,L)         . 

/^  h+2  h+3/2 

It  reme-ins  now  to  describe  the  formation  of  E.  and  of  dt 

J 

To  calculate  the  former  quantities  with  relations  (2.12)  we  need  first 
,j^b+3/2^  According  to  the  usual  stability  argnments  we  must  have 

^E    ^    ,  3rp.^/2      . 
At  y 

thia   is  equivalent  to 

/  Ob 

It  turns  out  that  the  right-hand  member  of  (3»36)  has  its  small- 
est value  at  the  shock.  Thus  v,  J^  >      d'&/  da     may  be  evaluated  there. 
This  yields  our  condition 
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iA.^'^   .  _^  _iB 


^  ^  ^  ^p     /  *  ^ 


r 


where  /   is  a  constant  <  1.  Our  formula  is  then  this: 

h+'^/2  '^  h+2 

Given  this  dt   '   we  now  proceed  to  the  calculation  of  B.    with  the 

help  of  (2.12).  To  this  end  it  is  necessary  first  to  calculate  E.  ,  i.e. 

J 

E       =  E(jds        ,  t  ). 

J 

This   is  accomplished  by  the  use  of  the  quadratic   interpolation  formula 
(3.19)   with 

f,  =^\  u  =  j(2  V  +  y^). 

J  J 

Then  (2.12)   yields  directly 

e/"2  (J    =1,   2,    ...,L-1)         . 

This  completes  the  induction.  E.    are  given  by  (3.31) >  (3-32);  p, 
by  (3.31^);  i/-'^  by  (3.27);  ^^^   by  (3.27);  D^-^^/^  by  (3.27);  D^'^/^ 
by  (3.2);  t^""^  by  (3. 27);  dt^-'^/^  by  (3.37)  and  dt^''^/^  ^y  (3.2). 

In  relation  (3.26)  we  adopted  what  may  at  first  seem  to  be  a 

^  ^  T^b+l/2 
quite  arbitrary  definition  for  when  two  approximants  to  D     are  near 

enough  to  each  other.  We  say  a  few  words  now  in  explanation  of  this 
formulation. 

Eelation  (3.9)  implies  that  if  (^  D  is  an  error  in  D,   <i"  (p  -  Pq)  , 
the  corresponding  error  in  p  -  p^,  is 


II- u 


and 

cf(p  -  Pq) 


2D 


P-P„  D^-    rv, 


iD. 


With  the  help  of  (3 '26)   thia  becomes 

cf  (P  -  Pq)  2D 


p  -  p^     ^  D.(yP0  )V^ 
At  the  'beginning  of  the  calculation  the  ratio 

is  about  7  and  it  approaches,  as  we  know,  assnaptotically  to  1.  In  the 
former  case  we  see  that 


(J(P  -  Pq) 


<  2£ 


and  in  the  latter 

»  -  % 

Thus  relation  (3.26)   guarantees  that  the  relative  error  in  the  shock 
over-pressure  determined  by  ,D  does  not  exceed  2£   for  large  D  and  £ 
for  D  near  (   ^ vJ  _fc)        '   ^^^   therefore  in  all  cases  the  relative 
error  in  the  determination  of  the  shock  over-pressure  is  about  2£  .  It 
is  important  to  notice  that  the  relations  (3.'+2)  and  (3.^3)  are  concerned 
with  over-pressure r:  since  they  are  the  relevant  quantities  in  the  compu- 
tation. As  will  be  seen  the  calculation  starts  with  this  over-pressure 
being  about  100  atmospheres  and  continues  until  it  is  about  0.02 
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atmospheres.  Evidently  it  is  the  relative  precision  of  this  quantity 
which  is  of  interest. 

We  close  this  section  vith  a  discussion  of  the  meaning  of  the 
iterative  relation  (3o25) .  As  we  saw,  for  each  given  D  there  are  deter- 
mined two  pressures  P/^+^/gW  P(h+l/2)  ^^^'^^  ^^®  equal  if  and  only  if 
the  D  is  the  correct  value.  Thus  we  desire  to  solve  iteratively  the 
functional  equation 

^(^)  =^[1x^1/2]  ^^^  -^(h+1/2)  ^^)  =° 
for  Do  As  is  well-known  an  error-aquaring  technique  for  solving  such 
relations  is  that  of  Wewton-Eaphson.  According  to  this  method  if  D  is 
an  approzimant  to  the  solution  of  (3«^^)j  then 

D^  =  D  _  _4^   (  ^'(D)  =  d^^/dD) 

is  a  better  one.   In  fact  if  D^  is  the  exact  solution,  i.e.  if 


then 
D 


^D^  =  ((D-D  •)  ^'(D)  ~  (  Z\(D)  -  Z\(D  )))/^'(D)~^^^  (D-D  )  = 
0        U  ^  2.£\ '  (D) 


Thus  the  error  (d  -  D  )  is  proportional  to  the  square  of  the  previous 
error  (D  -  D_) ,  the  factor  of  proportionality  being  A" /2.  A  '- 

The  exact  evaluation  of  ^^  '  is  quite  difficult  and  thus  instead 
of  using  precisely  (3^5)  we  have  used  (3-25).  In  that  formula  the 
esHPression 

(k^-k-.l^)/(k-l^  -^^^ 
is  introduced  as  an  approximation  to  -  l/  A   '  (jjP)  • 
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k"     The  Initial  valxiea » 

If  the  pressure  ratio  p/'p  at  the  shock  is  infinite,  it  is 
known  that  the  original  problem  (1.1)  -  (1-9)  has  an  exact  analytic  solu- 
tion. This  has  been  shown  by  the  senior  one  of  U3.   (Cf.  Erode,  Op^  Cit.) 
The  solution  is  a  similarity  one  in  the  sense  that  p  and  w  =  E/5^  are 
expressible  as  functions  of  the  variable 

Z  =  r/^   . 
It  is  not  difficult  to  show  that 
gi   =  a  t2/5„ 
To  express  p  and  B  it  is  convenient  to  introduce  a  parametric  representa- 
tion for  r,  E,  p. 


r=  ^^^^^^^/JJtil)    yJ^-^-^^;t^/2(r^/)  ]  /sr2-rjar^/) 


■-tl    r  -    'i/f^-lf-l-  12. 


R^gi^  f^n-xJ^/r^z]3r-J3fz-rjx^  r2f^/)  T  ^a-nrs^-o 


7-f 


With  the  help  of  the  equations  we  can  calculate  initial  values.   It  is 
known  that  for 

p/p  ~  100  atmos. 
the  equations  ('4-. 3)  yield  a  quit-?  good  description  of  the  situation. 

These  equations  are  somewhat  cumbersome  to  handle  in  a  machine 
calculation  so  we  have  replaced  them  with  the  following  equivalent  dif- 
ferential system.  In  this  system  r  has  been  replaced  by  s. 
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dK  7? 


/  -s^ 


ds  Zs:    [  f-^x 

dx      ^  x(/^x)i3(Z'f)X  -h  (Zf-^nJ 
ds  2s(f06^-^  Zyt  f)  > 


are  the  diffsrential  equations  and 
E(/)  =  ^, 
^(  ^)  =  1, 

are  the  initial  conditions.  The  integration  proceeds  from  a  =  ;/  to  s  =  0. 
5"  Preliminary  remarks  on  the  calculation. 

In  this  section  we  discuss  some  of  the  relevant  details  of 
the  calculation. 

As  we  mentioned  in  the  introductory  section  it  is  known  that  the 
analytical  results  (4.3)  give  a  close  approximation  to  the  actual  physi- 
cal situation  for  shocks  whose  pressure  ratios  are  in  the  neighborhood 
of  100  atmospheres.  Accordingly  it  was  decided  to  start  the  computation 
with  the  initial  shock  pressure-ratio  as  100  atmospheres  and  quite  arbi- 
trarily with  the  initial  shock  radius  as  l/2  unit;  also  arbitrarily  the 
pressure  and  density  of  the  undisturbed  gas  as  unity.  We  have  taken 
If  =   lA.  Hence 
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The  conditions  (5-1)  in  conjunction  with  {k.2)    determine  the  time  t  . 

With  these  initial  values  the  equations  C^-A)  ,  (iv-5)  were  inte- 
grated from  the  ahock  inward  to  yield  E,  p  at  the  time 
t°  =  0.0182575 


as  functions  of  jds  ,  where 


Thus 


d3°  =  1/800. 


K°,,  P°,  Wo     (j  =  1,  2,  ...,  L) 


y   ^  J-1/2 
are  fonned.  By  definition 

B%  =   0,  p°_-,^;2  =  p°^/2  . 

The  first  two  of  equations  (^V.if)  were  then  integrated  with 

where 

dt"^/^  =  dt-^/^  =  0.000086, 

D  =  2^/5t. 
In  order  to  tabulate  the  resulting  function  against  jds  the  equations 
were  first  integrated  from  J  ~^   to  ^7  °  -  ds  .  This  provides  starting 
values  for  E,  p  at  Lds*^.  From  this  point  inward  to  1  ds  the  calcula- 
tion proceeds  normally.  Thus  R  '  is  obtained  and  from  this  and  E  , 

J 

we  find 


J 

This  procedure  yields  the  starting  values  of 


T<iPr'iAl^U^-'U^'',^'''^"''^t'''''-''^^ 


they  arise  from  equations  (5-1-)  -  (5- 
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The  numerical  ■nrocedure  uaefl  to  carry  out  the  intearatior:  vas  the 
well-known  one  of  Runge-Kutta.  Three  different  intervals  of  integration 
were  used.     They  were 

-  As   =  1/32,000,   1/16,000,    1/8,000. 
We  compared  the  values   of  E  ,,   p      /^,  E"     as   obtained  by  tiie  tliree   inte- 
grations and  give  bslow  the  relative  errors.     Let  the  quantities  with 
"bars"   over  them  refer  to  those  obtained  with  the  interval  of  integra- 
tion    /As/   =  1/32,000;  with  "tildes"   over  them  the  interval    jA  s/    = 
1/16,000  and  with  no  symbol  the  interval    /A  3/    =  l/8,000.     Then 

^^  =  7.0/. I0-'  ^   M^  =  s.eZ.io-\  1^  =  ^,,7w.-^ 

The  data  corresponding  to     /A  s/    =  l/8,000  were  accepted  as  the 
initial  values  for  R,  p,  P- 

To  completely  characterize  the  problem  we  need  the  initial  data 
(5.1),  (5.3),  (5.9)  together  with  values  for  the  /^  of  (3«5?')  and  L. 
We  carried  out  three  different  calculations  determined  by 

pq  =  0,  L  =  100,  r=  1/2, 
p^  =  1,  L  =  100,  r=  1/2, 
pq  =  1,  L  -  80,  r=  2/5. 

The  results  of  these  caluclations  are  discussed  in  more  detail  in  what 
follows . 

The  actual  algorithm  can  best  be  described  inductively.   Suppose 
given 

fi.f.,ds)  di^-"^  di'^\  ^"'^^'>'^<<,^-lf '. 
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Tlie  modus  procedendi  is  then  this.  We  first  seek  the  value  D   '  .  This 
is  obtained  by  the  iterative  procedure  outlined  in  §3  above.  The  ap- 
proximant  ,D  is  given  by 

and  D  by  relation  (3.8).  The  iteration  then  proceeds  as  described  in 
§3;  the  principal  formula  being  (3.25)-  This  iteration  terminates  when 
the  relation  (3.26)  is  satisfied.  The  value 

vas  used  throughout  almost  all  calculations.  Near  the  end  of  the  prob- 
lem (5-11)  it  was  altered  to 

The  number  of  iterations  of  (3.25)  required  to  satisfy  (3-26)  with 
£  determined  by  (5-13)  was  about  three  independent  of  the  shock  over- 
pressure. To  test  "noise  level"  of  this  scheme,  i.e.  to  measure  the  ac- 
cumulation of  rounding  errors,  we  evaluated  the  D  corresponding  to  the 
shock  over-pressure  of  50,  starting  with  two  quite  dissimilar  sets  of 
initial  estimates.  The  relative  disagreement  was  about  one  part  in  a  mil- 
lion. For  smaller  over-press\ires  this  disagreement  was  even  smaller. 
Thus  our  choices  (5.13),  (5'13')  yield  results  above  the  "noiae  level". 
Having  calculated  D  "*"  '  ,  we  next  form 

and  then  with  the  help  of  (3 '37) 
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dt  ^^3/2^ 
It  is  interesting  to  notice  that  the  ratio 
dt^^3/2/,^h+l 

varies  only  by  a  small  amount  throughout  the  entire  calculation.  For 

an  over-pressure  of  about  ks   it  is  .Ok,   for  an  over-pressure  of  1.0^4- 

J  ) 

it  is  .08  and  for  an  over-pressure  of  .02  it  is  .03°  The  three  cor- 
responding values  of  ds    are 

.002,     .Qk,  2.92. 

These  quantities  having  been  determined,  we  are  enabled  to  cal- 
culate with  the  help  of  formulae  (2-9),  (2»15)  the 

_  h+i   _  -£.+1 

Then  with  the  help  of  the  intexpolaticn  formula  (3<>19)  we  form 

^  h+1  ^  h+1   A  h+2     . .   .   .       _ . 
j   '   ■^^-1/2'   1        \;^   -  L,,   X.,    > .. ,  L,) , 

and  this  completes  the  induction  since 

.h+1   .h   ,,h+l/2 
t    =  t  +  dt   '   . 

The  details  have  already  been  discussed  in  §3  above. 

There  is  one  important  deviation  from  this  program  we  must  de- 
scribe  at  this  time.  The  functions  E,  E  are  not  adeq\iately  represent- 
able  by  quadratic  approximations  near  s  =  0  to  permit  of  the  interpola- 
tion procedure  just  discussed.   In  the  table  below  we  give  the  values 
E_  ,  E,  ,  . .  . ,  B/-  to  show  this.  We  tried  to  run  the  problem  (5O10) 
as  just  described  and  found  that  quite  soon  it  reached  a  point  where 
E.  I  >  E,  "^  for  j  near  2.  This  physical  absurdity  stopped  the  cal- 
culation. 

To  rectify  this  situation  without  resoriiing  ~to   elaborate  methods 
we  recalled  that  E,  B  are  nearly  proportional  to  s  ' '  near  s  =  0,  as 
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can  be  seen  from  formulae  (^.3)  and  thus  that  R ' , '^'  are  nearly  linear. 
This  is  shown  in  the  following  table.  Accordingly^ we  first  raised  to 
the  power  7  the  functions  E,  R.  and  then  performed  the  interpolation 
for  j  =0,  1,  ...,  5.  For  j  >  5  the  functions  R,  R  are  sufficiently 
slowly  changing  to  permit  local  quadratic  approxinations. 


''o" 

0 

■'I" 

. 27^9^ 

«/ 

. 303^9 

< 

.32149 

\° 

.33^86 

< 

.3^^558 

.27^9^ 
.02855 
.01800 
.01337 
. 01072 


-.2I1559 
-.01055 
-.00^63 
-.00265 


i^.V 


0  7  -h  1»1877^I0'- 

(R^)  '  1.1877  :i-  10  -k.2   X  10  ' 

0  7  ^k  ^°^^-^5  3-  10"'' 

(E^)'  2.3712x10^  .    -5.3x10' 

0  7  -k  ^-1782  X  10" 

(E^  )  '  3.5^94  X  10  ,    -6o2  X  10  ' 

0  7  ^k  1°1'^20  X  10" 

(E,  )  '  1^.7211+  X  10  •  .        -6.9   X-  10  ' 

r,   7  ).  1.1651  X  10" 

(E  ") '    5-8865  X  10"^ 

In  a  typical  time  step  there  are  27'<-6  multiplications  and  1790 
divisions.  Since  the  only  variable  length  portion  of  the  code  is  the 
iteration  at  the  shock  and  since  the  number  of  these  iterations  is 
virtually  a  constant  this  count  of  the  multiplieationa  and  divisions 
is  quite  accurate.  Assuming  that  a  multiplication  reqviires  720  jus 
and  a  division  1100  lis  on  our  machine,  we  find  that  3«95  seconds  are 
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spent  in  each  time  step  on  multiplications  ar.d  dxTiaionso  An  actual 
"clocking"  of  the  time  spent  on  such  a  step  gives  10. IT*  Thus  the  solu- 
tion time  is  2.57  times  the  time  spent  in  multiplying  and  dividing. 

An  examination  of  the  solutions  [k.^)    shows  that  problem  (5»10) 
is  such  that 


^  h  /  ^-7v  h    h   /  h 


h 

are  independent  of  h,  i.e.  of  t  ,  To  teat  the  accuracy  of  the  code 
and  of  our  niimerical  procedure  we  accordingly  undertook  the  solution  of 
problem  (5. 10)  and  printed  out  the  data  (5»l6)<.  The  calculation  started 
with  a  shock-pressixre  of  100  atmospheres  and  stopped  with  a  shock-pressure 
of  about  hd.     To  measure  the  accuracy  of  the  results  notice  in  formulae 
{k.3)    or  (^.5)  that  the  shock  pressure  is 

r     25 (f-^ I)     i'        Z5fr*i)  -^ 

and  thus  that 

p  0(}  =   const. 
We  used  the  relative  change  in  this  quantity  p  "^  as  an  overall  figtore 

of  merit  for  our  code.  We  found  in  our  calculation  this  relative  error 

-k 
to  be  about  5.25  x  10  . 

The  results  R.'^,  p.'^Wp,  E ^■^■'-  (j  =  0,  1,  . . . ,  L)  were  printed 

out  periodically.  It  was  arranged  that  for 

p  >  2.0I+ 

this  took  place  whenever  a  12.5^  change  in  the  shock  pressure  occurred 

and  in  the  contrary  case  whenever  a  12.5^  change  in  the  shock  over- 

pressvire  occurred.  For  the  entire  problem  between  k^   and  i*-7  time  steps 

made  up  such  a  group. 
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The  problem  (5-11)  was  then  undertaken,  and  to  obtain  some  evalu- 
ation of  the  overall  accuracy  of  the  results  we  carried  out  (5. 12).  This 
yielded  an  estimation  for  the  overall  error  including  that  arising  from 
the  interpolation  process.  To  compare  the  results  of  the  two  calcula- 
tior^  we  used  the  results  of  the  problem  (5. 11)  to  obtain  the  values  of 
the  shock  over-pressures  at  the  shock  radii  obtained  in  the  problem 
(5.12)  and  then  compared  these  over-pressures  to  those  actually  obtained 
in  (5.12).  This  yielid  a  set  of  relative  errors  in  the  over-presaures 
which  increase  slowly  as  indicated  below. 

p°l  re lat ive  error 

50.00  1.8  X  lO'ji 

10.00  4.8  X  10"7 

5.00  8.2  X  10"7 

2.00  k.6   X  10"7 

1.00  5-9  X  10'^ 

.60  1.3  z  10":^ 

.30  1.8  X  10'^ 

.05  4.9  X  10"^ 

.02  7«9  X  10' 

The  comparisons  were  made  in  this  manner.  It  is  known  --  Of. 

Figure  I  and  Table  I  below  --  that  the  shock  over-pressure  and  shock 

radius  are  related  by  a  relation  of  the  form 

P  =  1  =  A  ^  -^ 
where  n  is  a  slowly  varying  function  of  y(    .  Further  as  we  indicated 
above  the  tabulations  we  obtained  from  the  machine  computation  occur 
at  12. 93^  changes  in  the  shock  over-pressure.  Thus^it  seems  reasanable 
to  perform  the  comparisons  discussed  above  logarithmically.  If  p  -  1, 

^   ;  p  -  1,  ^     are  two  consecutive  points  on  the  p  -  1,  ^  cui've 
aa  determined  by  the  solution  of  problem  (5'11)  and  p  -  1,  -A  is  a 
point  on  the  corresponding  cujrve  for  (5' 12)  with 
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then  our  interpolation  procedure  is  this.  Form 

log  (p  -  1)/(P2  -  1)  =  -  n  log  5^/  ^2-- 
where 

^  ^  .  log  (P3_  -  1)/(P2  -  1) 

Then 


where 


P     -/ 

and  we  have  accordingly  an  eatimate  for  the  relative  error  in  question. 

After  having  obtained  tabulations  of  the  functions 
E(r,  t)   p(r,  t) , 
we  required  tabulations  of  the  function 

p  =  P(E,  t). 
To  do  this  we  needed  to  invert  for  each  t  the  function 

r  =  E(r,  t^) 
and  substitute  the  resulting  function  of  R  into  p(r,  t  ).  This  was  ac- 
complished by  choosing  a  unit  dR  and  then  solving  the  equation 

for  u  as  a  function  of  k;  with  this  u  the  corresponding  p  can  be  found  aa 
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To  solve  (5»l8)  for  us  we  used  a  Nevton-Eaphson  type  formula.   If  u''  ia 
an  approximation  to  the  solution  of  (5.I8),  then 


,/  _  _L  /  V  ,     byU  ^  C 


where 


A  starting  approximant  is,  of  course 

Uq  =  -  c/Pb 

To  avoid  obtaining  too  many  values  P,  for  a  given  t  several  choices 

K  / 

were  made  for  the  unit  dR.   In  fact,  it  was  chosen  in  accordance  with 
the  relation 

d^  =  1/2  n  for  50/2''"^  >  ^  >  50/2"  . 

6.  The  results. 

We  have  arranged  the  results  of  our  calculations  in  a  number 
of  Tables  and  Figures  which  are  appended  ?iereto. 

In  Figure  I  we  have  displayed  the  shook  over-pressure  as  a  func- 
tion of  the  shock  radius.  Since  the  relation  between  these  quantities 
is  as  indicated  in  (5.17)  above  it  was  found  convenient  to  plot  the 
curve  on  log- log  paper. 

In  Figure  II  the  time  of  the  arrival  of  the  shock  at  a  given 
radius  is  plotted  against  that  radius.  This  plot  is  also  on  log-log 
paper. 

In  Table  I  followsing  Figure  II  we  have  displayed  the  shock 
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radius,  the  exponent  n  of  (5 '18),  the  dynamic  over-pressure,  and  the 
over-pressure  for  a  number  of  values  of  the  overpressure,  spaced  roughly 
1256  apart.   (The  dynamic  over-pressure  is  the  ordinary  one  plus  (I/2)  P   u  , 
i.e. 

P-l=p-l+  {1/2) f    u^ 
vhere  p,  ^  ,  u  are  the  pressure,  density  and  velocity  at  the  shock.)   In 
Table  II  the  shock  velocity  and  time  of  arrival  are  displajred  at  the  same 
over-pressure  as  in  Table  I. 

Figures  III-VII,  inclusive,  are  graphs  of  the  over-pressure  as 
a  function  of  Eulerian  radius  at  various  typical  times. 

Figures  YIII-X7I,  inclusive,  are  graphs  of  the  over-pressure  as 
a  function  of  the  time  for  various  values  of  the  radius. 

Figure  XVII  is  a  graph  of  the  length  of  the  positive  phase  as  a 
function  of  the  shock  radius. 


Figure  I   -  11-28. 
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Shock  Badlua 

EzpotMnt 

Orerjireiaur* 

OfTpnaBwr^ 

.5004  74 

0 

.  00000000 

1  07 

.267  0  67 

100 

.0064  54 

.525  760 

2 

.97610288 

9  2 

.673  4  12 

86 

. 36  16  48 

.5  50  567 

2 

.9  5953521 

8  0 

.8  91937 

75 

.34  64  92 

.576  675 

2 

.954189  71 

7  0 

.583003 

65 

.708  4  40 

.604  178 

2 

.94841455 

6  1 

.56  1980 

57 

.27  4879 

.6  33179 

2 

.9  41785  90 

5  3 

.668  4  92 

49 

.89  5875 

.66  3  795 

2 

9  33830  94 

4  6 

.762  940 

43 

.44  09  28 

696  158 

2 

.9  2521041 

4  0 

.72  1328 

37 

.79  4133 

.7  304  16 

2 

.  915001  53 

35 

.436805 

32 

.85  56  19 

.765  9  44 

2 

903955  20 

30 

.906  5  12 

28 

. 62  2687 

.803  632 

2 

89157209 

26 

.9  33  221 

24 

.9110  26 

8  4  3  688 

2 

.8  77  259  37 

23 

.449499 

21 

.6576  36 

8  86  3  50 

2 

.86122505 

2  0 

.395646 

18 

.806757 

.9  31894 

2 

.  84  357076 

17 

.719069 

16 

.309  279 

.9  80  637 

2 

.82361186 

1  5 

.374122 

14 

12  25  94 

0329  50 

2 

.801359  90 

13 

.32  0  561 

12 

209  147 

089  2  61 

2 

.77529243 

1  1 

.52  390  0 

10 

53  6  768 

150070 

2 

74  6299  86 

9 

.953026 

9 

0764  39 

215  961 

2 

7  174  50  71 

8 

.579  0  99 

7 

.80  1273 

287  616 

2 

68311671 

7 

.3796  39 

6 

.  69  0302 

.365  8  31 

2 

6  4  5955  62 

6 

.333  4  63 

5 

.  72  37  53 

.451 543 

2 

61593589 

5 

.4  18  756 

4 

88  1288 

5  52  2  70 

2 

.44926905 

4 

.6  12605 

4 

14 1609 

.6  61705 

2 

4  78565  83 

3 

.908110 

3 

498  135 

.7  88  86  3 

2 

40611748 

3 

.282131 

2 

92  9  4  44 

.9  37  649 

2 

.33963248 

2 

.728914 

2 

429  9  90 

082  148 

2 

.27  4766  82 

2 

.32  00  60 

2 

.  06  32  34 

24  3  4  33 

2 

.21204  4  74 

1 

.968011 

1 

.749  343 

428095 

2 

.14611345 

1 

.6  59  866 

1 

.4762  18 

,6  35  562 

2 

07982242 

1 

.397  344 

1 

.24  4782 

8  68  984 

2 

014  12061 

1 

.174  533 

1 

.049  2  11 

137  9  43 

94882107 

0 

.982414 

0 

88  10  91 

448  715 

8  8  4  200  52 

0 

.818071 

0 

737  4  72 

8  08  693 

.82050117 

0 

.6  78613 

0 

61 5528 

112  652 

767708  57 

0 

.589  463 

0 

.537405 

447  3  41 

7  244  32  12 

0 

.512301 

0 

.469  579 

8  15  8  05 

683085  33 

0 

.445  591 

0 

.410702 

22  1359 

64406084 

0 

.3  87  949 

0 

.359  5  80 

6  80  28  8 

6  06828  92 

0 

.3  3690  5 

0 

31  40  56 

.186  450 

5714  33  15 

0 

.2  92  975 

0 

27  46  33 

.744  546 

.5  38  38082 

0 

.2  55  145 

0 

24  04  62 

,376  744 

.5  07  102  78 

0 

22  1747 

0 

.2100  88 

075  204 

4  77  630  12 

0 

.193  042 

0 

. 18  3801 

8  67  4  77 

44988167 

0 

167  748 

0 

. 16  04  78 

744  576 

4  23766  93 

0 

146  0  26 

0 

.14  03  11 

7  15  3  62 

3  99  560  94 

0 

127  338 

0 

.12  28  48 

8  17  898 

376805  19 

0 

110  856 

0 

107  3  52 

04  08  18 

3554  10  03 

0 

096  675 

0 

09  39  39 

397041 

335449  04 

0 

084  4  48 

0 

08  23  12 

9  39  538 

316914  70 

0 

073  646 

0 

07 1987 

6  54  0  20 

2  99  194  85 

0  . 

06  4  328 

0  . 

06  30  39 

5  59  4  49 

2825244  6 

0  . 

0  5t>  ii77 

0  . 

05  5274 

7  300  74 

2669918  3 

0  . 

04  9  150 

0  . 

04  8  374 

148  072 

25238925 

0  . 

04  2  988 

0  . 

04  2387 

8  414  34 

23844940 

0  . 

037  650 

0. 

037184 

915013 

225092  74 

0 

0  32  919 

0  . 

03  2559 

346  797 

211886  80 

0  . 

028  822 

0. 

028  5  44 

178  348 

200299  34 

0  . 

02  5  2  66 

0. 

02  50  51 

456  0  60 

18893540 

0  . 

0  2  2  17  5 

0 

. 02  2008 

346  4  89 

170925  21 

0  . 

019443 

0  . 

019  314 

8  19  7  20 

151133  84 

0  . 

017  086 

0  . 

016985 

Table   II   -  11-31. 


Ttne  of  Shock 

Overpreaaure 

Shock  Velocity 

Arrival 

100 

006  4  54 

0 

.11018518 

0  . 

0  18  300 

66 

361646 

0 

10248609 

0  . 

02  0  6  82 

75 

346  4  9  2 

0 

.09582056 

0  . 

0  2  3  186 

65 

7  08  4  40 

0 

06958244 

0  . 

0  2  6  0  0  :> 

57 

27  4  879 

0 

.08  3743  57 

0  . 

0  2  9  182 

49 

8  95  87  5 

0 

.07827838 

0  . 

0  32  766 

4  3 

4  4  0  9  28 

0 

07316359 

0  . 

0  36  814 

37 

794  133 

0 

06837614 

0  . 

0  4  13  9  1 

32 

8  55619 

0 

06389581 

0  . 

0  4  6  5  77 

28  . 

6  2  2  687 

0 

05978898 

0  . 

0  5  2  3  28 

24  . 

9  110  26 

0 

05594035 

0  . 

0  5  8  848 

21  . 

6  57  6  36 

0 

05233466 

0  . 

0  6  6  2  55 

18  . 

8  06  757 

0 

04895723 

0  . 

0  7  4  6  8  8 

16  . 

309  2  7  9 

0 

045794  25 

0  . 

0  8  4  312 

14 

122  594 

0 

04283353 

0 

0  95  3  24 

1  2 

209  147 

0 

04006367 

0 

.107  959 

10  . 

5  36  76  8 

0 

0374  75  4  9 

0 

.12  2  501 

9  . 

076  4  39 

0 

0  35059  56 

0 

13  9  2  88 

7  . 

8  0  12  7  3 

0 

03280477 

0 

158  7  29 

6  . 

6  90302 

0 

0  3070  5  64 

0 

18  13  2  1 

5  . 

72  3  753 

0 

02875501 

0 

20  7  6  60 

4  . 

8  812  8  8 

0 

02693983 

0 

2  38  4  74 

4  . 

14  16  09 

0 

02523  8  72 

0 

2  77  168 

3  . 

498  135 

0 

02365959 

0 

3  2  2  0  0  1 

2 

9  2  9  4  44 

0 

02217055 

0 

3  77  5  76 

2 

429990 

0 

02077496 

0 

4  4  6  978 

2 

06  3  2  34 

0 

01968726 

0 

5  18  4  88 

1 

74  9  34  3 

0 

01870618 

0 

6  0  2  6  0  2 

1 

4  76218 

0 

01780860 

0 

70  3  86  1 

1 

24  4  78  2 

0 

01701099 

0 

8  2  3  158 

1 

049211 

0  . 

01630661 

0 

.963420 

0 

8  810  91 

0 

01567581 

1 

.13  1777 

0  . 

7  37472 

0 

01511610 

1 

.333817 

0  . 

6  15  5  28 

0 

01462407 

1 

.576  110 

0  , 

5  37  4  05 

0 

01429995 

1 

.786380 

0  . 

46  9  57  9 

0 

01401247 

2. 

0  2  2  9  0  1 

0  . 

4  10  7  02 

0 

0137580  6 

2 

2  88  3  60 

0  . 

3  59  580 

0 

0  135  33  28 

2 

5  8  5  6  6  1 

0  . 

314  0  56 

0 

01332992 

2 

9  27  4  44 

0  . 

274  633 

0 

01315127 

3 

30  9  8  28 

0  . 

240462 

0 

01299444 

3 

7  36  84  2 

0 

,210088 

0 

01285343 

4 

22  6  119 

0 

.18  3  8  01 

0 

012730  12 

4 

7  72  249 

0 

.160478 

0 

.  01261972 

5 

397  433 

0 

140  311 

0 

.  01252347 

6 

0  95  226 

0 

.12  2  8  48 

0 

.0124  39  53 

6 

8  73  112 

0 

107  35  2 

0 

.01236455 

7 

76  2  2  19 

0 

.093939 

0 

01229930 

8 

.731101 

0 

.082  312 

0 

.  01224244 

9 

.836444 

0 

.07198  7 

0 

. 01219174 

1  1 

0  9  9  125 

0 

.063039 

0 

.01214762 

1  2 

508  0  4  5 

0 

.055  2  74 

0 

.  01210921 

14 

0  7  9  1  9. 3 

0 

.048374 

0 

.  012074  97 

1  5 

8  74  377 

0 

.042387 

0 

. 0120451P 

17 

8  79  4  40 

0 

.037  184 

0 

.  01201924 

20 

118  000 

0 

.032  559 

0 

.  01199613 

2  2 

6  77  7  80 

0 

.02  8  544 

0 

.  011976  03 

2  5 

54  10  22 

0 

.025051 

0 

. 011958  51 

28 

7  4  2  810 

0 

.022008 

0 

.  011943  24 

32 

32  2  317 

0 

019  3  14 

0 

.  011929  70 

36 

4  19  4  66 

0 

016  9  8  5 

0 

. 011917  98 

4  1 

00  9  7  19 

Figure  III   -   11-32. 
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Figure  V  -  II-3^. 
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Figure  X  -  11-39- 
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Figure  XIII   -   11-^2. 
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Figure  XVIII   -  Il-k-J. 
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B.  CONTINUED  FRACTION  EXPANSIONS  OF  ALGEBRAIC  NIMBESS* 
0.  Summary. 

Every  real  irrational  number  x  has  a  unique  infinite 
simple  continued  fraction  expansion 


=  l\'  ^2'  ^-y   °°°] 


V 

where  the  a.  are  integers,  positive  for  i  >  1;  conversely  every  such 
infinite  simple  continued  fraction  represents  a  unique  real  number.   (A 
rational  number  has  two  trivially  different  finite  simple  continued 
fraction  expans ion.s  <. ) 

Among  the  irrational  numbers  the  quadratic  surds  are  distinguished 
as  those  whose  simple  continued  fraction  expansioi^s  are  eventually  peri- 
odic --  for  example,  Z^'^   --=  [1^,2,1,2,...]   =  [l,l?2  ]  . 

The  question  naturally  arises  as  to  whether  there  are  any  evident 
peculiarities  in  the  simple  continued  fraction  expansions  of  other 
algebraic  numbers.  As  material  for  such  an  inquiry,  extensive  expan- 
sions of  gi%'en  algebraic  numbers  were  desiredo  Codes  for  producing  them 
have  been  written  a.nd  used  on  the  IAS  computer,  and  work  is  in  progress 
on  additional  codes. 

The  first  code  written  dealt  with  the  case  next  beyond  the  quad- 
ratic, that  is,  with  cubic  irrationalities.  With  modifications  of  this 

1/3 
code,  the  continued  fraction  expansions  of  2  '   to  2232  terms  and  of 


*  A  portion  of  this  problem  has  been  submitted  in  a.   preliminary  note 
to  MATEEMATICAL  TABLES  AND  OTHER  AIDS  TO  COMPUTATION  for  publication. 
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3  '   to  500  terins  wers  obtained,   (The  iLodifications  chiefly  inrolved 
methods  of  checking  the  computations  in  the  face  of  ciemory  difficulties 
related  to  the  high  niemory  reference  repetition  rates  characteristic  of 
this  code..) 

A  code  is  nov  being  written  which  generalizes,  from  3  to  n,  the 
degz'ee  D  of  the  polynomial  defining  the  algebraic  irrationality,  and 
which  also  incorporates  error-detecting-and-correctiug  facilities o  Addi- 
tional production  of  results  is  planned  when  this  more  general  code  is 
available . 

The  results  so  far  obtained  have  been  examined,  without  any  obvi- 
ous peculiarities  being  no-ciced.   In  addition,  certain  statistics  have 

been  tabulated,  such  as  the  frequency  distribution,  sum,  and  geometric 

1  /3 
mean  of  the  first  2000  partial  quotients  of  2  ''  . 

Khint chine  and  L^yj  have  given  some  probabilistic  results  on  the 

distributions  of  the  sum  S  (x1  and  of  the  geometric  mean  G  (x)  of  the 

n'  ■  ^  n^  ' 

first  n  partial  quotients,  and  of  the  nth  root  of  the  denominator  of  the 
nth  convergent,  of  random  numbers  x  uniformly  distributed  between  0  and 
lo  The  results  for  the  sura  are  of  an  essentially  weaker  type  than  those 
for  the  two  other  variables  mentioned. 

For  example,  for  almost  all  x  --  i.e.,  except  on  an  exceptional 
set  of  measuxe  zero  --  G  (x)  approaches  an  absolute  constant  K.  One 
may  ask  whether  or  not  a  given  number  belongs  to  this  exceptional  set 
(evidently  most  if  not  all  quadratic  surds  do  belong  to  it) .  However, 
the  theorems  do  not  in  themselves  provide  numerical  estimates  of  the  dis- 
tribution functions  of  G  (x) ,  so  that  the  interpretation  of  a  G  (x) ,  com- 
puted from  a  finite  number  of  partial  quotients  of  a  given  z,  is  still 
open. 


I 
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1.  The   jorLtiaued  fraction  algorithm  for  real  algebraic  irra- 
tionalitiea. 

It  is  desired  to  find  the  simple  continued  fraction  expan- 
sion of  a  real  irrational  x  defined  as  a  root  y  =  x  of  an  integral  poly- 
nomial P(y)  of  degree  D,  We  will  suppose  that  x  is  >  1  and  is  the 
largest  root  of  P(y)  (if  not,  the  first  step  of  the  algorithm  would  need 
a  small  modification) . 

Suppose  that  the  property,  "there  is  given  an  integral  polynomial 
P  (y  )  ,  of  degree  D  in  the  indeterminate  y  ,  having  y  =  x  =  [_a  ,a   ,.o.j 
as  its  sole  root  >  1/'  holds  for  n  =  h.  Then  the  integral  part  of  x,  is 
the  desired  h-th  partial  quotient,  a,  ,  characterized  as  the  largest  inte- 
ger a  for  which  sgn  F(a)  ^   sgn  F(+  oo).  Let  7x^  =  a^  +  ■'"/^h+l*  ^^®"^  ^^ 

is  easily  shown,  (since  x,  =  a.  +  l/x.  ,)  that  if 
^      h    n    '  n+1' 

^h+i^yh+i^  ^  ^Li  •  \K  ^  ^'\^i^ 

then  the  above  property  holds  for  n  =  h  +  1.  Letting  PTiy-,)  be  P(y) ,  all 

the  a  and  P  are  then  defined  inductively-  The  above  procedure  is  the 
n     n  J  f 

algcrithm  used  for  determining  the  partial  quotients. 

2.  Outline  of  the  basic  code. 

Each  step  of  the  algorithm  is  seen  to  consist  of  three 
t3rpes  of  operations  with  polsmomiala  F  of  degree  D; 

(1)  determination  of  t:he  largest  a  (=  a,  )  for  which  sgn  F.  (a)  /  sgn  F.  (+  od  ) 

(2)  translation  of  the  independent  variable  by  a.  , 

(3)  inversion  of  the  independent  variable. 

The  compiiter  mist   be  coded  to  accomplish  a  sequence  of  these  steps. 

As  the  induction  index  h  increases,  the  coefficients  of  the  poly- 
nomials tend  to  increase,  and  quickly  exceed  the  40-binary-digit  word 


i 
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length  of  the  computer.  Hence  it  is  necessary  to  define  a  representa- 
tion of  large  numbers  as  p~tuples  of  computer  words.  At  any  stage  of 
the  computations  for  an  irrationality  of  degree  D,  we  need  D+l  such 
multiple-word  length  numbers,  and  it  may  be  useful,  but  is  not  essential, 
to  have  an  additional  one  as  a  working  space.  It  is  convenient  to  repre- 
sent these  D+l  or  D+2   numbers  with  a  common  multiplicity  pi  to  avoid 
waste  of  machine  time  in  the  early  part  of  the  compu'^ations  (since  com- 
putation time  is  nearly  proportional  to  p) ,  p  is  made  variable,  starting 
at  1  and  increased  only  as  necessary,  as  the  coefficients  grow.   (Even- 
tually the  coefficients  grow  to  occupy  all  the  storage  available  for 
them.  This  limits  the  lengt.h  to  which  any  particular  non-periodic  ex- 
pansion can  be  computed,  since  the  computer  has  a  finite  memory.) 

A  chief  requisite  for  the  continued  fraction  code  is  the  provision 
of  subroutines  for  doing  the  necessary  arithmetic  on  multiple  word-length 
(MWL)  numbers o  Since  MWL  arithmetic  demands  codewise  lengthy  methods  of 
accomplishing  mathematically  simple  operations,  it  pays  in  the  present 
codes  to  keep  the  b8,sic  multiple-word-length-arithmetic  subroutines  (MWL 
SEs)  free  of  multiplication  orders  (for  time  economy) ,  and  limited  in 
variety  and  complexity  of  purpose  (for  coding  and  memory-space  economy), 
relying  on  controlling  subroutines  to  call  the  MWL  SEs  as  needed  to 
synthesize  less  simple  operations. 

The  operations  chosen  for  realization  as  MWL  SEs  ares   left  shift 
and  right  shift  (i.e.,  multiplication  by  specified  positive  or  negative 
—  or  zero  --  powers  of  2) ,  and  addition.  Each  MWL  SE  is  entered  with  a 
three -parameter  word  specifying  the  MWL  operands  —  or  operand  and  num- 
ber of  shifts  --  and  the  location  of  the  MWL  result. 

The  determiriation  of  a.  (^  1)  is  effected  binary-digit-wise,  by 
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ascertaining  first  the  position  of  the  most  sig-nificant  digit  (1)  of  x,  , 
and  then  the  presence  or  absence  of  1  in  each  successively  less  significant 
position  (down  through  the  units  position) ;  or  stated  differently,  by  suc- 
cessively determining  and  subtracting  from  z,  the  highest  remaining  power 
of  2  (down  through  2  )•  The  translation  by  a,  ,  from  P,  (x,  )  as  a  polynomial 
in  X,  to  P..(ai^  +  \)   as  a  polsmomial  in  z,  (=  x,  -  a,  )  is  made  stepwise, 
by  powers  of  2,  as  the  corresponding  digits  of  a,  are  determined » 

Each  test  of  a  polynomial  F(u)  to  determine  whether  sgn  F(2  )  / 

k 
sgn  F(+  oo),  followed  by  a  translation  by  2  in  case  the  inequality  holds, 

is  in  practice  carried  out  by  testing  a  polynomial  F*(v)  =  F(2  v)  to 

determine  whether  sgn  F*(l)  /  sgn  F*(+  oo),  followed  by  a  translation  by 

1  in  case  the  inequality  holds.  Starting  from  a  P,  ,  the  exponent  k  of 

the  scale  factor  2  is  first  increased  (from  0)  until  the  highest  power 

of  2  in  a_  is  found,  then  decreased  to  0  while  the  remaining  powers  of  2 

are  found. 

The  inversion,  converting  Pj^Ca-u  +  z,  )  into  P>,.T(yh.T)  = 
7-^^1   '  ^■^(\   ■•■  ■'-/y'h+l) '  wtiere  z  =  l/yv,^i'  consists  merely  in  interchang- 
ing the  coefficients  of  each  two  terms  whose  degrees  in  z  add  to  D. 

In  this  manner  the  inductive  step,  of  determining  a,  and  Fj^  -, 

from  P,  ,  is  reduced  to  a  combination  of  the  elementary  poljmomial  opera- 
h 

tiona: 

(a)  replacing  F(u)  by  F*(v)  =  F(2v) 

(b)  replacing  F(u)  by  F*(v)  =  F(v/2) 

(c)  replacing  F(u)  by  F*(v)  =  F(l+v) 

(d)  replacing  F(u)  by  F*(v)  =  v   •  F(l/v) 

(e)  testing  whether  F(l)  ^  F(+  w   )• 

These  operations  on  polynomials  can  in  turn  be  broken  down  into  operations. 
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by  the  three  MWL  SEs,  on  their  coefficients.   (In  the  original  code 
these  breakdowns  were  made  in  extenso  for  the  case  D  =  3-   In  the  more 
general  code  in  progress,  and  in  an  intermediate  check-code  already 
written  and  used,  they  are  done  inductively,  and  are  valid  for  arbitrary 
D.) 

Thus  the  inducti-^re  step  of  the  algorithm  is  constructed.  With 
its  incorporation  into  a  repetitive  process,  and  the  addition  of  input, 
output,  and  other  peripheral  facilities  (such  as  the  periodic  read-out 
of  the  coefficients,  from  which  the  code  could  be  re-started  at  a  later 
time),  the  basic  program  is  complete. 

3-  Checking  facilities. 

For  results  to  be  reliable  it  is  essential  to  have  methods 
of  checking  them.  At  first  with  a  new  code  the  concern  is  with  the  cor- 
rectness of  the  coding;  then  in  production  the  concern  is  with  the  ac- 
curacy of  the  machine  calculations.  Among  the  possible  methods  of  check- 
ing are  the  following: 

(1)  Duplication;  that  is,  making  supposedly  the  same  computation 
twice,  and  examining  whether  the  results  are  identical; 

(2)  Using  two  somewhat  different  methods  which  should  give  the 
same  results,  and  examining  whether  the  results  are  identical; 

(3)  Using  a  redundancy  check;  that  is,  some  property  (not  used 
directly  in  performing  the  codes)  which  should  hold  among  the  intermediate 
or  final  results. 

In  the  continued  fraction  codes,  a  large  fraction  of  the  computa- 
tion time  is  spent  in  the  inner  loops  of  the  ^WL  SEs,  where  corresponding 
operations  are  performed  on  one  segment  after  another  of  the  two  or  three 
MWL  numbers  involved.  These  loops  have  few  orders,  with  no  multiplications 
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or  divisions,  and  therefore  have  a  high  memory-reference  repetition  rate. 
A  well-known  hazard  ("read -around")  of  cathode-ray  tube  memories  is  that 
frequent  reference  to  the  same  memory  locations  may  endanger  the  contents 
of  nearby  locations  which  are  not  concurrently  referred  to;  one  way  of 
conventionally  avoiding  this  danger  is  to  isolate  the  high-reference -rate 
locations,  by  a  buffer  zone  of  unused  locations. 

In  spite  of  such  isolation  of  the  inner  loops  of  the  MWL  SEs,  the 
continued  fraction  code  has  at  times  been  more  than  usually  subject  to  ma- 
chine errors,  some  of  which  have  even  been  duplicated.  The  evidence 
strongly  suggests  that  the  errors  occur  in  the  memory,  and  are  related 
to  the  high  memory-reference  repetition  rate,  but  their  nature  remains 
obscure;  they  are  not  of  the  type  usually  designated  as  read-around,  and 
have  not  been  so  troublesome  with  other  codes. 

In  the  continued  fraction  code  (more  than  in  some  other  codes)  a 
certain  type  of  machine  error  is  propagated  with  increasing  effect  through 
all  later  calculations;  so  that  no  errors  can  be  tolerated. 

As  a  check  on  machine  operation,  method  (1),  while  wasteful  of 
time,  is  simple,  and  is  an  effective  detector  of  certain  machine  errors. 
No  results  have  been  accepted  without  at  least  this  check.  Against  the 
above-mentioned  errors,  which,  repeated  in  identical  circumstances  (even 
though  the  overall  error  rate  was  low) ,  duplication  offers  no  protection. 
It  is  no  check  at  all  on  the  coding. 

Method  (;:;)  is  also  wasteful  of  time;  but  by  imposing  a  different 
pattern  of  computations  or  of  digits,  or  a  restated  logic,  it  could  show 
up  some  types  of  coding  errors,  and  in  the  detection  of  machine  errors 
is  safer  than  duplication.   In  this  problem  it  was  early  applied,  pri- 
marily as  a  code  check,  by  slightly  changing  the  basic  code  so  that 
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instead  of  the  polynomials  P  ,  there  were  produced  in  one  case  3  •  P^ , 

n  h. 

in  another  case  k  •  P  .  The  agreement  or  proper  relationship  between 
results  was  a  partial  confirma.tion  of  the  correctness  of  the  MWL  SEs 
(which  were  also  checked  by  cases) ,  since  compensating  errors  would  be 
unlikely  to  occur  in  such  different  numbers . 

A  method  of  T.jrpe  (3)  ,  which  is  appropriate  for  the  continued  frac- 
tion codes,  is  to  check  the  arithmetic  on  the  MWL  numbers  by  performing 

analogous  arithmetic  on  their  residues  to  some  modulus,  which  for  con- 

,  39  , 
venience  has  been  taken  as  (2  -1) ,   and  comparing  the  results.  This 

process  is  fast,  compared  to  the  main  computation,  since  it  deals  with 
single -word -length,  instead  of  WfTL,    coefficients.  Although  it  does  not 
check  all  phases  of  a  computation,  it  does  give  positive  protection 
against  the  greatest  hazard,  the  change  of  a  single  binary  digit  in 
some  coefficient  (it  could  also  show  up  some  types  of  coding  errors). 
A  check  code  using  this  process  (for  general  degree  D)  was  written,  and 
used  to  check  the  results  already  produced.  The  generalized  main  code 
now  being  written  incorporates  the  check,  concurrent  with  the  main  com- 
putation, so  that  if  an  error  is  detected,  the  computation  may  be  re- 
started immediately  from  a  just  previous  point;  this  will  be  more  ef- 
ficient than  a  post-mortem  method. 

Additional  results  will  be  produced  when  this  more  general  and 
efficient  code  is  available. 
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C.   MAGNETOHTDEODINAMIC  THEORY  OF  SOIAE  SPICm'ES 

Thia  problem  is  being  compu+,ed  for  Dr.  Schlliter,  Dr.  Schwarzs child, 
and  Mr.  Lautman  of  the  Princeton  University  Observatory. 

Spicules  or  jets  are  relatively  small  (6"  of  arc)  disturbances 
of  the  solar  chromosphere  which  are  characterized  by  an  increasing  rate 
of  outward  propagation  and  a  short  lifetime  of  the  order  of  two  minutes. 
Knowledge  is  urgently  needed  on  the  nature  of  spicules  since  it  is  now 
thought  that  they  transport  energy  to  the  corona. 

In  this  theory  the  spicule  begins  as  a  pressure  variation  at  the 
base  of  t;he  chromosphere.  In  the  absence  of  a  magnetic  field  this  would 
create  an  ordinary  spherical  sound  wave  in  the  Isothermal  chromosphere. 
However,  the  effect  of  a  magnetic  field  on  a  disturbance  in  the  highly 
conductive  gas  gives  rise  to  the  phenomenon  known  as  magnetohydrodynamic 
waves  (see  Alfven,  Cosmical  Electrod3mamic3,  Oxford,  1950).  The  velocity 
of  propagation  of  these  waves  is  inversely  proportional  to  the  square 
root  of  the  density,  hence  in  the  very  high  density  gradient  of  the 
chromosphere  the  propagation  velocity  of  the  wave  will  increase  greatly 
as  it  progresses. 

The  basic  equations  of  this  problem  3-re; 
1.  E  =  -vxH     (l=(j*(E+vxH)  with  <3^  =  oo 

X  E   = 

3.        ^  z  H  =  i^  TT   J  (-— ^  «  J) 


k. 


dv  T 

nJl^jxH-VP-  (^^  <P  (equation  of  motion) 
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5-  -^:4- —  =  -  "V  .  (  P  v)     (equation  of  contimiity) 

The  aasumptions  made  ares 

1.  Plane -parallel  chromosphere. 

2.  Barometric  density  gradient  (  p    ~  e~^  ). 

3.  Small  vertical  extent  (  'v' CD   =  g  =  const.)- 

h.     Axial  symmetry  (H>,  =v>,  =0,  E  =E  =J  =J  =0). 
'      T        z     s     z     s 

5.  Homogeneous  vertical  magnetic  field  of  value  H„. 

6.  All  changes  adiahatic. 

Assumption  6.  leads  to  another  basic  equation,  the  adiabatic 
equation; 

at  dt 

Under  these  assumptions  the  basic  equations  are  expanded  and 

linearized,  and  the  following  non-dimensiorial  variables  defined: 

p 
r=-^s,   J  =  ^  z,  •J;'=-'fc,   "="B"=  const. 

C  C  C  Jo 

2                        H^ 
W  =  ^^v^,   V^^^v^,   B=-| =  B^e^ 

Sg  ^»  C  ^„ 

P  =  indisturbed  pressure 
O     =  indisturbed  density 
The  linearized  equations  then  take  the  form: 

7-  j^  Br:|i".i-si].if^(Bw) -i^(i) 
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^  ^r  ^ 


2 

10.      _  Z  =  _s;)    £_   p    +  (-0   -  1)  V 


Equations  7-10  may  be   combined  to  yield  the  final  differential 
equations   of  the  problem; 


U.      ^MB.V, 


■^r 

,    (BW)    4 

+  (-J 

-)^ 

br 

"^^). 

r  >  ^  .  2W 

t/  '-  c)y  d)r  or  ()  j  dr 

The  right-hand  sides  of  (11)  and  (12)  are  converted  to  finite- 
difference  relations  by  the  equations  r  =E^r  (0<E<  15)  and  y  = 

Ya  y  (O  <  Y  <  17)   by  means   of  the  usual  definitions   of  the  space- 

Y  ^  y    ^  y 
derivatives.  B  becomes  Be     .  e     is  computed  at  the  beginning 

of  the  problem,  and  with  each  succeeding  vertical  increment  is  multiplied 

A   y  ^  y 

by  e       "^   to  give  its  proper  value-     The   computation  of  e       "'is   simplified 

by  means   of  the  following  expansion 

e  ^  y  =   ((((^^  +  1)-^  +  1)^^  +  1)    ^  y  +  1)        . 
1+3  2 

The  initial  conditions  are: 
_  ^   /  2^  2. 

P  =  P  e~   ^   ^  '  (pressure  "hump") 

W  =  V  =  0  (all  velocities  equal  zero) 

Therefore  initially  7-  becomes 

'^  W  ^   1^  "^  P 
'^^     gr  "^r 

and  8.  with  the  aid  of  9*  becomes 


i 
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1)  V 


1  "^P 


g  c>y 
We  are  now  able  to  compute  the  initial  values  of  W  and  V  at  all 
points  in  the  field.  To  simplify  the  boundary  conditions  the  wave  ia 
put  in  a  "box",  i.-e.  horizontal  velocities  vanish  at  vertical  boundaries 
and  vertical  velocities  vanish  at  horizontal  boundaries.  The  singu- 
larities in  11.  at  r  =  0  are  removed  by  means  of  an  expansion  valid  in 
this  region. 

W  and  V  for  each  time  step  are  computed  from  the  formulae 

•T+I/2   'T-l/2  ,   .  _  'Jr    ,  ^T+1   ^T    ,  ~  'T+l/2   ^ .     .   ^^ 
X   '   =  X   '   +  at:  X  and  X    =  X  +  ^  ^  X   '  .  Since  m  the 

*T=0 
final  equations  the  only  time -dependent  term  is  X  and  X    is  known  from 

the  initial  conditions,  it  is  a  very  simple  matter  to  compute  each  suc- 
ceeding time  step  from  the  one  before. 

The  32  X  32  place  memory  of  the  computer  is  divided  for  the  pur- 
pose of  this  problem  into  three  fields  as  shown  in  the  following  diagram. 


17 

Order  Stage 

V  and  V 

double 

stored 

W  and  W 

double 

stored 

0     r     15      r     31 

At  each  point  in  the  field  At  W  and  AT  V  are  computed  and 

*T     '-T  T+l     *T+1 

added  to  W  and  v  to  give  W    and  V   ,  these  quantities  being  stored 

in  the  second  half  of  each  word-  After  all  points  have  been  computed  we 

T  'T+1 

have  W  in  the  first  phase  and  W    in  the  second  phase.  The  same  is 

true  for  V.  The  computer  is  then  made  to  sweep  through  the  whole 
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bottom  part  of  the  memory,  multiplying  the  second  pahae  of  each  word 
by  AZ^  and  adding  the  reault  to  the  first  phase.  This  then  gives  us 
the  new  W  and  V  at  all  points.  This  two-step  procedare  ia  more  efficient 
than  a  single  scan  over  the  field  giving  an  immediate  computation  of  the 
new  W's  and  V's,  the  second  method  requiring  complicated  shift  procedures. 

Since  there  is  not  enough  room  in  the  memory  (the  magnetic  drum 
was  not  in  operation  at  the  time  this  problem  was  coded) ,   a  converaion 
code  was  not  incorporated,  the  results  being  punched  out  in  binary  form. 
The  results  are  then  converted  in  a  separate  operation.  Five  quantities 
are  required  at  each  point,  W,  W,  V,  V,  and  P   ,  however,  we  are  only 
able  to  store  the  first  four  for  all  points.  It  would  be  very  time- 
consximing  to  place  the  five  quantities  in  temporary  storage  and  punch 
out  point -by-point  as  the  computation  proceeds.  Therefore,  the  compu- 
tation  of  o  by  9*  is  incorporated  in  the  conversion  code.  Thus  it  is 
only  necessary  to  punch  out  \^   cards  containing  II52  quantities  as 
against  288  cards  containing  ihhO)   quantities.  Punchouts  are  made  every 
twelfth  time  step.  In  addition,  a  supplementary  one-card-punchout  is 
made  containing  the  volume  integrals  of  P   for  each  of  the  12  time 
steps.  This  is  a  check  to  see  that  the  enclosed  mass  remains  constant. 

Two  integrations  have  been  performed  and  others  are  pending.  One 
integration  was  performed  with  B  =  0  (no  magnetic  field)  and  the  expected 
spherical  wave  resulted.  Another  integration  was  performed  with  a  mag- 
netic field,  but  the  results  have  not  been  analyzed  as  yet.  The  pending 
integrations  will  test  the  coarseness  of  the  time  increments  for  the 
second  integration  above  as  well  as  observe  the  effect  of  different  mag- 
netic field  strengths  both  on  the  above  density  gradient  field  and  on  a 
homogeneous  density  field. 


II-61. 


D.   CAICULATION  OF  TEE  EMERGY  BAND  STRUCTOBE  OF  IRON  WITH  APPLICATION 
TO  THE  THEOEY  OF  FEREOMAONETISM. 

The  fundamental  problem  of  solid  state  phsrsics  ia  the  calculation 
of  the  energy  values  and  wave  functions  of  electrons  in  a  solid.  On  this 
basis  it  should  be,  in  principle,  possible  to  explain  all  the  observed 
electronic  behavior  of  solids.  But  because  a  solid  contains  so  very- 
many  electrons,  a  rigorous  calculation,  which  would  have  to  consider 

23 
some  10   mutually  interacting  particles  for  a  sample  of  ordinary  size, 

is  utterly  out  of  the  question. 

The  most  common  approximation  which  is  made  in  order  to  reduce 
the  problem  to  tractable  dimensions  is  that  the  wave  function  for  the 
solid  can  be  written  as  some  simple  combination  of  "one -electron"  func- 
tions --  functions  which  depend  only  on  the  co-ordinates  and  the  spin 
of  one  particle.  This  wave  function  has  to  satisfy  the  requirements  of 
anti-symmetry  imposed  by  the  Pauli  exclusion  principle;  the  wave  func- 
tion of  the  system  must  be  anti-symmetric  under  the  interchange  of  the 
coordinates  and  spin  of  any  two  electrons.  A  form  of  the  wave  function 
which  satisfies  this  requirement  is  the  determican'^al;  the  wave  function 
is  written  as  a  determinant  of  one-electron  functions,  one  for  each 
energy  states 


fn. 


U,a,,5,)    .  .  .  M,(Xo,5„) 


U„(Z,,^,)     ...       ^ri(^rij>S^ 


The  quantities  x....x  are  the  coordinates  of  the  n  electrons,  and  the 
^  1    n 

quantities  s..  ...s  are  their  spin. 


^ 
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According  to  the  minimum  theorem  of  quantum  mechanics  the  best 
possible  wave  function  of  the  system  will  be  obtained  by  minimiziag  the 
expectation  value  of  the  energy  with  respect  to  small  variations  of  the 
one-electron  function  u.(x.,  s.)  (providing  these  functions  remain  nor- 
malized and  orthogonal) .  This  variation  procedure  leads  to  an  equation 
for  the  one-electron  functions: 

The  Hamiltonian  operator  can  be  written  as 

in  which  i  and  j  refer  to  the  position  of  the  electrons,  and  <f  refers 
to  the  atoms.  V^(  r" .)    is  the  potential  energy  for  electron  i  in  the 
field  of  atom  ', nucleus  plus  core,  (if  there  are  low  lying  states  which 
are  not  altered  by  presence  of  other  atoms).  The  equations  obtained  by 
the  minimization  are  the  Hartree-Fock  equations 

H(x^)  contains  the  kinetic  energy  for  electron  ^  plus  its  potential 
energy  in  the  field  of  all  the  nuclei.  The  quantity  (  -fc.)  is  the  one- 
electron  energy.  The  last  term  on  the  left,  or  exchange  term  vanishes 
unless  the  spin  i  =  spin  j. 

It  has  been  shown  that  the  Hartree-Fock  equations  can  be  regarded 
as  ordinary  Schroedinger  equations  for  the  motion  of  electrons,  each 
electron  moving  in  a  slightly  different  potential  field,  which  is  com- 
puted by  electrostatics  from  all  the  changes  of  the  system,  positive  and 
negative,  corrected  by  the  removal  of  an  exchange  charge,  equal  in  mag- 
nitude to  one  electron,  surrounding  the  electron  whose  motion  is  being 
investigated.   In  a  solid,  the  periodic  arrangement  of  the  atoms  causes. 
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byBloch's  theorem.a  like  periodicity  of  the  one-electron  functions  which 
are  solutions  of  the  Hartree-Fock  equations.  To  each  wave  function  cor- 
responds an  energy  level,  and  these  levels  are  grouped  in  the  energy 
hands . 

The  Hartree-Fock  equations  can  be  used  as  a  basis  for  the  theory 
of  Ferromagiietism.  Consider  an  unmagnetized  crystal  in  which  we  reverse 
the  spin  of  one  electron,  from  -  to  +.  Since  there  are  now  more  pairs 
of  electrons  of  positive  spin  than  of  negative  spin,  the  exchange  term, 
which  reduces  the  coulomb  repulsion  of  the  electrons,  tends  to  lower  the 
energy  of  the  crystal.  On  the  other  hand,  it  is  necessary  to  place  the 
electron  whose  spin  is  changed  into  a  higher  energy  state,  for  all  the 
lower  states  of  that  spin  are  filled.  Ferroraagnetism  will  occur  when 
the  decrease  in  exchange  energy  is  greater  than  the  increase  caused  by 
the  operation  of  the  exclusion  principle:  the  Fermi  energy. 

In  principle,  the  Hartree-Fock  equations  should  be  solved  by  the 
method  of  the  self -consistent  field.  This  has  so  far  proved  impossible 
for  the  following  reason;  assuming  a  starting  potential,  the  Hartree- 
Fock  equations  yield  a  set  of  wave  functions.  It  is  necessary  to  com- 
pute a  new  potential  and  compare  and  average  it  with  the  first,  and  re- 
peat the  process  until  the  potentials  from  the  n-1  and  the  n-th  stages 
agree  with  an  allowed  error.  But  the  calculation  of  the  second  potential 
requires  a  summation  over  all  the  occupied  states,  and  our  methods  of  solv- 
ing the  differential  equations  in  a  solid  are  not  sufficiently  good  so  as 
to  give  accurate  solutions  except  for  states  of  special  sjnmnetry. 

One  has  to  use  his  physical  intuition  as  to  what  corresponds  to 
a  reasonable  crystal  potential,  including  the  effect  of  exchange  terms, 
and  try  to  solve  the  differential  equation  for  as  many  states  as  possible 
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using  this  potential.  One  coramon  way  to  obtain  a  starting  potential  ia 
to  superpose  the  charge  densities  from  Hartree  fields  for  free  atoms 
placed  at  the  proper  points  in  the  crystal  lattice o  An  approximate 
exchange  potential  can  be  obtained  in  the  same  way  using  Slater's  free 
electron  approximation  in  which  the  exchange  potential  is  proportional 
to  the  cube  root  of  the  charge  density.   If  it  is  desired  to  examine  the 
possibility  of  Ferromagnetism  for  the  system  it  is  necessary  to  prepare 
exchange  potentials  for  the  magnetized  and  unmagnetized  states » 

Given  a  starting  potential,  one  way  of  solving  the  Hartree-Fock 
equation  is  through  the  orthogonalized  plane  wave  method.  This  approach 
expands  the  wave  function  in  terms  of  a  complete  set  of  fimctions  periodic 
in  the  crystal;  plane  waves  which  have  been  made  orthogonal  by  the  Schmidt 
process  to  the  (assumed  known)  low  lying  states  of  the  atomic  coreo  An 
orthogonalized  plane  wave  is  defined  as: 

where  J/q   is  the  volume  of  a  unit  cell  of  the  crystal,   /7^  runs  over 

all  atomic  sites,  ^j(l^~l^t)    is  "t^®  J-th  core  function  on  the  site 
/ 
r   and 

The  wave  function  u,  which  is  classified  as  belonging  to  the  i~th  ir- 
reducible representation  of  the  wave  vector  \,   is  written  as: 
./ 


k 
We  must  form 


u 
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where  H  is  the  appropriate  Hamillonian  from  tne  (simplified) Hartree-Fock 
equations.  The  ajniibol  (  )  repreasata  scalar  product.  From  (7)  we  ob- 
tain a  set  of  linear  equations  for  the  coefficients  a^,  which  must  lead 
to  a  sscular  equation  of  the  form: 

det  [(X, ,  /V/J  -  €^(K,  X,)]=0 

We  will  include,  of  course,  only  &   finite  number  of  terms  in  the  series 
{h) ,   hoping  that  the  series  will  essentially  have  converged.  The  process 
of  solving  this  secular  equation  is  greatly  simplified  by  group  theoretical 
considers-tions  according  to  the  ssnmnetry  of  the  function  u^. 

This  work  has  5.nvolV9d  a  calculation  of  the  energy  values  and  wave 
functions  for  the  body-centered  form  of  iron,  according  to  the  method 
whose  principles  hs-ve  been  sketched.  The  greatest  part  of  the  work  con- 
sists in  the  calculation  of  the  matrix  elements  (  X. ,       hj '^l     )  ^^-^ 
{     X        y,  )}   a'^'3  the  solution  of  the  equation  (8)  for  l6  different 
states  of  high  sjinmetry.  Equation  (8)  is  equivalent  to  a  problem  of 
simultaneous  diagonalization  of  two  Hermitian  matrices,  and  it  is  in  this 
form  that  the  solutions  have  been  found.  Since  it  is  desired  to  investi- 
gate the  ferromagnetism  of  iron,  exchange  potentials  have  been  prepared 
for  three  different  spin  arrangements,  corresponding  to  a  magnetized  and 
an  unmagnetized  state.  The  solution  of  (8)  for  the  l6  states  has  to  be 
repeated  for  each  exchange  potential.  This  problem  would  be  completely 
impossible  without  the  use  of  a  high  speed  computer. 

After  the  energy  values  are  available  for  certain  states  of  high 
symmetry,  it  is  necessary  to  interpolate  to  find  the  energies  of  more 
general  states  according  to  some  procedure  and  then  to  compute  a  density 
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of  states-curve  (number  of  states  per  unit  energy  as  a  function  of 
energy)  for  each  exchange  potential.  The  comparison  of  the  density  of 
states-curves  should  reveal  the  tendency  toward  Ferromagnetism,  in 
that  a  decrease  in  the  energy  of  the  system  should  result  on  passing 
to  the  ferromagnetiG  from  the  unmagnetized  state. 

We  give  in  the  following  paragraphs  a  brief  description  of  the 
mathematical  m'ittod  used  to  effect  the  calculations  carried  out  on  this 
problem.   It  is  clear  that  one  desires  to  find  the  proper  values  and  the 
corresponding  functions  of  the  matrix  A  -  /I  B  where  A  is  symmetric  and  B 
is  positive  Hermitian.  In  a  prior  report  (Final  Report  on  Contract  No. 
riA-36-03^-ORD-1023)  we  discussed  a  method  for  the  case  in  which  B  is  the 
identity  matrix.  From  the  mathematical  side  our  interest  in  the  present 
problem  lies  in  examining  the  situation  where  B  is  a  quite  general  Hermitian 
matrix.  To  handle  the  present  problem  we  replace  the  matrix  A  by  a  new 
aymmetric  matrix  C  such  that  the  proper  values  of  C  -  /\  I  are  the  same  ss 
these  for  A  -  /\  B  and  we  also  establish  a  procedure  for  finding  the  cor- 
responding proper  functions.  We  proceed  in  this  fashion. 

By  means  of  our  procedure  for  handling  s^Tnmetric  matrices  we 
express  B  in  the  form 

B  =  U*DU 
where  U  is  a  unitary  matrix,  D  is  a  diagonal  one,  and  U*  is  the  transpose 
of  U  (we  assume  all  quantities  to  be  real;  this  is  not  a  fundamental  re- 
striction; everything  that  is  being  done  is  still  valid  if  one  replaces 
symmetric  by  Hermitian  and  makes  the  obvious  emendations  in  the  text) . 
We  now  extract  a  reciprocal  square  root  of  B  in  this  manner 
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and  express  a  matrix  C  In  the  following  fashion 

C  =  D'^'^ATJ*!?-'^   =  (U*D'-^'^)*A(U*D'-'-'^). 
Consider  now  a  proper  value  A  and  a  corresponding  proper  func- 
tion z.  By  definition 

(A  -/\  B)   z  =  0. 
In  vhat  follows  it  is  convenient  to  define  a  new  matrix  V  for 

which 

-1/2 
V  =  U*D  ' 

and  a  vector  j  for  which 

y  =  v""''  X  =  D^'^   X 

Then  we  have 

(C  -  A  I)y  =  (V*AV  -  ?\l)j  =   V*(A  -  A  {Y*'^)Y'^)Yj  = 

=  V*(A  -  /\  JJ*li^'%^'^)Yy  =  V*(A  -  AB)Vy  = 

=  V*(A  -  Ab)  X  =  0. 
Thus  every  proper  value  /\  of  the  matrix  A  -  Ab  is  also  a  proper  value 
for  the  matrix  C  -  /\  I  and  the  corresponding  proper  vectors  are  related 
as  indicated  in  the  relation  (#) .  The  converse,  of  course,  is  equally 
clear. 

Thus  the  procedure  for  handling  the  matrix  A  -  Ab  has  been  re- 
duced to  that  of  handling  the  matrix  C  -  A  I.  First  the  matrix  B  is 
diagonalized,  then  C  is.  The  procedure  indicates  a  quite  convenient 
method  for  inverting  a  matrix.  Although  it  involves  somewhat  more  work 
than  the  Gaussian  elimination  method,  it  is  somewhat  more  stable  than  that 
method  and  indicates  precisely  where  the  loss  of  precision  occurs.  The 
inversion  of  a  unitary  matrix,  of  course,  in  volves  no  loss  of  precision. 
The  place  at  which  the  precision  is  lost  is  precisely  in  the  inverting 


*om 
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of  the  elements  cf  the  diagonal  matrix  D.  Furthermore,  the  loss  of  pre- 
cision is  encountered  in  a  quite  explicit  fashion  in  that  it  arises  frc 
the  shifting  which  must  he  performed  to  keep  all  elements  properly  scaled. 
Thus  one  obtains  the  inverse  matrix  and  an  exact  determination  of  its 
accuracy. 

Let  W  be  the  unitary  matrix  associated  with  C;  i.e.,  let 

w*cw  =  yi 

where  /I  is  a  diagonal  matrix,  and  let 

Y  =  W 
where  V  is  defined  above.  Then  if  our  results  are  absolutely  accurate, 
we  will  have 

AY  -  yA   =  0. 
Accordingly,  we  can  measure  the  error  in  our  calculations  by  evaluating 
the  size  of  AY  -  YA       relative  to  the  size  of  Y. 


PABT  III   -  METEOROLOGY 
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A.  Computations  with  a  three-level  model. 

During  the  year  several  further  tests  were  made  of  the  three- 
level  model*  using  both  the  R-  -  and  (T  -coordinate  systems.  These 
12  and  2^  hour  forecasts  were  for  cyclonic  developments  which  occurred 
in  November  1952  (two  cases),  February  1953,  and  November  1953-  A  fur- 
ther case  of  cyclone  filling  during  March  195^  is  presently  under  con- 
sideration. 

The  results  of  these  tests  were  considered  very  satisfactory. 
Figures  1,  2,   and  3  show  the  900  mb  results  for  the  forecasts  of  November 
24/25,  1952.  The  forecasts  were  made  using  the  0~  -coordinate  system 
with  initial  data  at  hOQ,   TOO  and  900  mb**.  From  these  figures  it  will 
be  seen  that  while  by  no  means  perfect  the  forecasts  are  as  good,  if 
not  better,  than  these  which  could  be  obtained  by  the  more  conventional 
methods  of  forecasting. 


*  For  a  description  of  this  model  see  Report  on  Contract  No.  N-T- 

ONE-388,  T.  0.  I  and  Final  Report  on  Contract  No.  DA-36-034-0RD-1023. 

**  The  initial  data  for  this  case  were  supplied  by  Dr.  Sverre  Petterssen 
of  the  University  of  Chicago  who  requested  the  forecasts  for  use  in  a 
study  of  cyclogenesis  currently  being  undertaken  at  Chicago. 
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B.  Experiments  vith  a  two -level  forecasting  model. 

DiLring  the  year  1952-53  some  forecasts  were  made  with  a  two-level 
model  using  data  at  the  300  and  TOO  mb  surfaces.*  These  forecasts,  while 
improving  on  the  barotropic,  or  one-level,  forecasts,  were  not  unduly  suc- 
cessful and  the  model  was  soon  superseded  by  a  three-level  model.  Subse- 
quently suggestions  were  made  by  other  investigators  that  if  data  had  been 
used  at  other  levels,  viz.  the  500  and  1000  mb  surfaces,  much  more  accurate 
forecasts  would  have  been  obtained.  As  a  result  of  these  suggestions  a 

2.  2. 

model  using  p-     (or  '^     in  tlse  notation  used  in  the  description  of  the 
three-level  model*)  as  vertical  coordinate  was  coded  for  the  machine.  The 
levels  at  which  the  forecasts  were  made  were  500  and  866  mb  but  the  initial 
866  mb  data  were  obtained  by  linear  interpolation  with  respect  to  pressure 
between  500  and  1000  mb  and  a  corresponding  1000  rob  forecast  obtained  by 
similar  extrapolation  from  the  forecast  500  and  866  mb  heights. 

The  equations  used  were  exactly  analogous  to  those  used  in  the  cor- 
responding three-level   6~  -system  forecasts.  They  involved  a  time  ex- 
trapolation with  the  resulting  need  to  compute  a  Jacobian  and  the  solution 
of  a  pair  of  simultaneous  linear  differential  equations  of  the  form 


=r  O 

=    o 


where  ^  ,  ^  are  functions  of  the  horizontal  position  coordinates  (x, 
y)  and  ^  and  ^  are  the  unknowns.  Solutions  were  obtained  for  a  19  x 
19  grid  using  p   =  constant  as  the  boundary  condition.  The  solutions 


*  See  report  on  Contract  No.  N-7-onr-388,  T.  0.  1,  and  Final  Eeport  on 
Contract  No.  DA-36-03^-0RD-1023. 
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were  obtained  157  the  same  iterative  procedure  used  for  the  three-level 
model.  Approximately  lyjo   "over -relaxation"  was  used. 

The  results  of  the  forecasts  made  with  data  from  November  1950 
and  November  1953  showed  no  significant  improvement  over  the  previous 
two-level  forecasts  and  the  model  was  therefore  not  used  for  further 
forecasts. 
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C.  Experiments  with  m.ilti -layer  non- linear  fore caa ting  models . 

Several  non-linear  models  were  experimented  with  diiring  the  year. 
The  accuracy  of  the  forecasts  was  unfortunately  affected  ty  an  incon- 
sistency in  the  phjrsical  assumptions  made.  This  inconsistency,  resulting 
from  the  introduction  of  the  geostrophic  assumption,  is  apparently  very 
difficult  to  remove  if  we  wish  to  retain  the  non-linearity.  Despite 
this  difficulty  the  computations  were  valuable  firstly  because  they 
brought  to  light  this  physical  inconsistency  and  secondly  because  they 
demonstrated  the  feasibility  of  solving  certain  non-linear  partial  dif- 
ferential eqviations.   In  order  to  demonstrate  the  method  of  solving  the 
non-linear  equations  a  brief  mathematical  description  of  the  non-linear 
version  of  the  tb^ee-level  model  is  given  below. 

In  the  report  for  last  3raar*  if;  was  shown,  how  the  equations  for  an 
n-lasrer  model  could  be  derived.  These  equations  are 


f 


=r   O 


(>4  -  !^A^> 


i^.  -A)'"''^ 


=    o 


■^^^. 


7^-/ 


I  [iX  \  %  (^--  -  <J"""/  =  ^ 


^- 


where    '^  >/     ~     ^  t/L  ll^-A  ~~  ^-^-n  J      ^"^^  ^^^  general  notation  is  as   in 

the  above  mentioned  report,  namelyj     <p    is  the  geopotential,  ^  Z.     ,  3>/I)-fc. 

the  individual  derivative  in  an  isobaric  surface,  ^  =  -^  -h    -^  V  ^ 


*  Ibid, 
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the  vertical  component  of  absolute  vorticity,  0^    the  U.  S.  Standard 
Atmosphere  potential  temperature  at  the  level  k.  The  pressure  interval 
p  =  0  to  p  =  P-.  is  assumed  divided  into  n  equal  subintervals,  ^/i-   , 
and  quantities  at  the  mid  points  of  each  interval  are  denoted  by  the  sub- 
script k  (k  =  1,  2,  . . . ,  n)  and  points  at  the  upper  and  lower  end-points 
of  each  interval  "by  k  -  1/2  and  k  +  I/2,  respectively. 

We  now  define  a  q,  by 
^k  ' 


and  therefore  the  n- layer  model  equations  (1)  become 


with  corresponding  equations  for  the  other  two  levels,  k  =  1  and  n. 

The  forecasting  system  then  is  very  similar  to  that  described  for 
the  linear  versioii,  the  only  difference  being  in  the  definition  of  q,  and 
thus  in  the  method  of  finding  the  p^    corresponding  to  a  given  field  of 

"^jL  .  The  computations  were  performed  for  the  case  k  =  3  using  the  0"  - 
coordinate  system.   In  this  case  the  three  aquations  to  be  solved  become 


^i  =  I3  (^-  -  '^') 


r 

These  are   clearly  a  very  non-linear  set  of  equations  for  the  ^   's. 

Having  indicated  the  actual  eqtiations  solved  it  is  convenient 
for  descriptive  purposes  to  now  ret^orn  to  the  general  form 
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Define 

Kow  assume  that  a  small  change  G_^    in  'p^   causes  a  change  o  ^/^ 
in  /t^ 

We  now  neglect  second  and  higher  powers  ot    &  f^ 


/3^    _  _£kl  _  ±^ 


Now  if  we  wish  to  red^Jice  a  "residual",  /C-^  ,  to  zero  by  only  altering 
'p^      we  must  choose  E-^     such  that  o^   =  ~^ik     y   i'©*  we  must  choose 
^11^  to  be 

Thus  if  <^M_       and  p^       are  successive  values  of  ^/i    during  the  it- 
erative process 
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Similar  eq.uations  clearly  hold  for  k  =  1,  n.  In  the  case  considered  in 
the  computations  the  tlaree   q\is.tions  were  solved  by  the  simultaneous  ap- 
plication of  procedure  (k)   for  all  three  values  of  k  at  a  given  point  (i, 
j)  the  points  being  scanned  in  the  (i,  j) -plane.   It  was  found  that  it 
was  not  sufficient  to  give,  as  in  the  linear  case,  the  criterion 

Mr'-  >^r/  <  ^ 

as  that  for  convergence.  Rather  it  was  found  necessary  to  stipulate  a 

certain  minimum  number  of  iterations  (six  in  the  cases  tested)  and  then 

to  require  the  usual  criterion.  No  attempt  was  made  to  use  "over-relaxation". 

The  applicability  of  this  method  clearly  depends  on  the  obtaining  of 
a  good  first  guess,  m ^     ,  in  the  iterative  sequence.  It  is  apparently 
sufficiently  good  in  the  meteorological  case  to  take  as  a  first  guess  the 
value  of  52)   obtained  by  linear  extrapolation  from  the  corresponding  values 
half  hour  and  one  hour  previously.  It  was  further  found  that  it  was  necessary 
to  restrict  ^^   to  positive  values.  If  any  negative  ones  appeared  in  the 
course  of  the  timewise  application  of  equation  (1),  these  values  were  re- 
placed by  zero's.  From  the  meteorological  view  point  this  assumption  was 
acceptable  as  the  validity  of  the  geostrophic  assumption  is  in  grave  doubt 
for  negatives  -^  's. 
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D.  Barotropic  forecasts  for  low  latitudes. 

In  conjunction  with  Professor  Eiehl  of  the  University  of  Chicago 
a  series  of  twenty  'barotropic  forecasts  were  made  for  areas  extending 
to  within  10°  latitude  of  the  equator.  Separate  200  and  50O  mb  fore- 
casts were  made  for  two  sumraer  series  and  5OO  mb  forecasts  for  a  spring 
series. 

The  only  code  differences  between  these  forecasts  and  the  baro- 
tropic forecasts  previously  made  arose  from  the  revised  formulae  neces- 
sary for  the  computation  of  t,he  sine  of  the  latitude  and  the  mapping 
factor  for  regions  of  low  latitude.  The  grostrophic  assumption  was, 
however,  still  ms.de  even  for  low  latitudes. 

The  forecasts  at  the  500-inb  level  consistently  brought  out  the 
major  circulation  featuji'ea  aad  in  many  cases  smaller-scale  features  were 
also  accurately  predicted.  Forecasts  at  the  200-mb  level  were  much 
poorer  and  showned  systematic  errors  which  resulted  in  some  very  marked 
failures. 

It  is  expected  that  these  predictions  will  be  of  considerable 
assistance  in  the  forfflijilation  of  more  accurate  models  for  numerical 
prediction  in  low  latitudes. 
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E.  Computation  of  vertical  velocitiea. 

In  the  forimilation  of  the  three-level  model  we  have  at  one  stage 
equations  involving  w,  the  individual  derivative  of    jl^    ,  aX   two  levels. 
In  the  normal  forecasting  system  it  is  convenient  to  eliminate  these 
two  variables  and  this  ia  in  fact  done.  Since,  however,  w  is  proportional 
to  the  vertical  velocity  and  is  therefore  highly  correlated  with  rain- 
fall amounts  it  is  interesting  to  know  its  distribution  at  various  stages 
in  the  forecast. 

If  the  levels  for  which  forecasts  are  made  are  numbered  1,  2,  and 
3,  w  can  be  evaluated  at  levels  1  l/2  and  2  l/2  by  use  of  the  relation 
ships 


where  c/,!/^   and  c/  , y   are  constants,  the  operator  ^r  =  -y^  -hY-V 
and  the  p   's  are  the  geopotentials  at  the  various  pressure  surfaces 
indicated  by  the  subscripts. 

In  order  to  obtain  the  time  derivatives  involved  it  was  necessary 
to  know  the  ^  values  at  more  than  one  time.   In  practice  it  was  found 
that  smooth  fields  of  w  were  obtained  if  the  time  derivatives  were  re- 
placed by  three  hour  centered  finite  differences.  Shorter  time  inter- 
vals resulted  in  too  much  round-off  error  and  therefore  ragged  fields 
of  w. 

Used  in  conjunction  with  a  three-layer  model  for  the  400,  TOO 
and  900  mb  surfaces  this  code  gave  vertical  velocities  at  560  and  805 
mb.  Exact  verifications  against  observation  is  impossible  because  of 
our  inability  to  measure  vertical  velocity  directly  by  observational 
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means.  In  the  cases  considered,  however,  the  results  agreed  with  those 
expected  from  other  considerations  and  also  agreed  fairly  well  with  ob- 
served rainfall  areas. 
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F.  Time  series  code. 

At  the  request  of  Dr.  Hans  Panofsky  of  Pennsylvania  State  College 
a  code  was  written  to  obtain  certain  auto-correlations  and  cross-correlations 
for  given  sets  of  observations.  The  observations  used  were  obtained  in  the 
course  of  an  investigation  of  small  scale  turbulence  in  the  atmosphere. 
The  code  is,  however,  very  general  and  can  be  used  to  treat  many  other 
types  of  time -series. 

The  input  to  the  code  consists  of  two  sets  of  initial  data,  --<^<.  and 
yi^;^  ,  i  =  1,  2,  ...,  N.  This  data  are  first  "pre-whitened"  yielding  sets 
of  data  Xl  i     /L      defined  in  the  following  manner: 

where 

for  i  =  n*,  n*  +  1,  ...  where  n*  <  ?• 

The  auto-correlation  and  cross-correlation  terms  are  then  computed 
for  -^  =  0,  1,  . . . ,  m,  0  <  m  <  128.  These  terms  are  defined  in  the  fol- 
lowing ways 
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I 


^('^  ^G^zTT^)  1  y.  X,.  -,^±^t1x1  ^. 


^TTt*'  t.--n*-t- 


The  code  then  goea  on  to  compute,  for  k  =  1,  2.  . . . ,  m  -  1,  the 
following  quantities: 


where 


and,  if  only  one  "pre -whitening"  parameter,  Q  f     ,  is  used. 
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For  k  =  0  and  k  =  m,  ^<1  ,  y     ,    Z4    ,  VJ     are  given  by  similar 
expressions  with  the  additional  factor  I/2  included. 

Several  sets  of  data  were  treated  in  the  way  outlined  above. 
These  sets  had  various  values  of  W  and  m  but  in  all  cases  to  date  O ^ 
has  been  taken  as  -  0-75' 
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G.  Objective  analyalg. 

All  methods  developed  to  date  for  numerical  weather  prediction 
require  input  data  at  a  regularly  spaced  set  of  grid  points.  Meteorological 
observations,  however,  are  made  at  an  almost  random  set  of  observing  sta- 
tions. We  thus  have  the  problem  of  obtaining  from  this  set  of  observations 
the  most  probable  values  at  the  grid  points  as  required  by  the  forecasting 
code. 

The  usual  method  of  solving  this  problem  is  to  plot  the  actual  ob- 
servations on  a  niap,  draw  a  series  of  isopleths  of  the  quantity  concerned 
at  suitable  intervals,  and  then  using  judicious  interpolation  to  read  off 
the  values  at  the  grid  points.  The  isopleths  are  drawn  using  continuity 
in  time  as  well  as  in  the  three  space  coordinates.  This  method  apparently 
depends  considerably  on  the  skill  of  the  analyst  who  draws  the  isopleths 
and  is  therefore  referred  to  as  "subjective"  analysis.  The  main  drawback 
to  this  method  is  not  its  subjective  nature  so  much  as  the  fact  that  it 
is  very  time  consuming. 

In  the  past  the  numerical  forecasting  models  have  been  tested  only 
on  very  carefully  aaalyTied  data  the  preparation  of  which  took  a  time  in- 
compatible with  any  routine  forecasting  system.  Two  questions  thus  arise; 
Is  so  much  care  necessary?  Does  the  subjective  analyst  do  anything  that 
a  computing  machine  could  not  also  do  in  a  considerably  less  time?  The 
only  way  to  find  answers  to  these  questions  is  to  propose  alternatives  to 
the  methods  used  bjr  the  subjective  analyst  and  then  to  test  these  alterna- 
tives by  using  them  in  the  preparation  of  forecasts.   In  this  section  one 
such  alternative  is  described  and  results  given  to  show  its  effects,  firstly, 
on  the  type  of  analysis  produced,  and  secondly,  on  numerical  forecasts  from 
data  so  prepared. 
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The  method  is  purely  two-dimensional  and  is  based  on  the  fitting 
of  a  second  order  polynoaiial  to  the  observations  in  the  area  inmediately 
surrounding  the  grid  point  at  which  the  value  of  the  pressure  or  height 
is  required.   In  this  description  we  will  be  concerned  with  the  height  of 
given  pressure  surfaces  expressed  in  terms  of  height  deviation,  D,  from 
standard  atmosphere.  The  polynomial  can  be  written 

Thus,  if  the  origin  of  the  coordinate  system  is  taken  to  be  the  grid  point 
under  consideration,  the  "D"  value  we  seek  will  be  <3.q^  .   In  order  to 
determine  the  ais  coefficients  Clij  we  clearly  require  six  independent 
observations  of  D  at  known  positions  (x,  y) . 

The  observations  give  us  wind  values  as  well  as  D  values.  These 
wind  values  are  incorporated  in  the  scheme  by  using  the  geostrophic  as- 
sumption. Thus,  one  wind  observation  gives  us  values  of  'h~J>/hs^       and 
'b'\)  /'du       at  a  point.  Provided  that  there  is  at  least  one  "D"  value 
within  the  region  being  considered  each  wind  observation  is  equivalent  to 
two  "D"  values  or  pieces  of  information.  The  errors  introduced  by  the 
geostrophic  assumption  are  not  random.  However,  due  to  our  lack  of  knowl- 
edge of  the  geostrophic  deviations  we  have  to  assume  that  the  errors  from 
this  souirce  are  random. 

Since,  in  general,  the  observations  from  nearby  stations  are  not  in 
exact  agreement  we  usually  require  more  than  the  minimum  number  of  pieces 
of  information  and  then  fit  the  second  order  polynomial  to  these  observa- 
tions by  the  m.ethod  of  Least  Squares.   Thus  we  choose  the  ^  ij   such  that 
we  minimize 
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where  i>^     ,   'ol>y^  /dzx.     ,   ^7)-k  /c)  'a        are  from  the  observations  and  the 
corresponding  primed  quantities  a-re  those  computed  from  the  polynomial 
(1).    ^^'2>     ^^   "^^^  standard  deviation  for  the  errors  in  the  -f~    observa- 
tions of  D  and  "^^  is  the  standard  deviation  of  the  S  wind  observations 
in  the  area  over  which  the  polynomial  is  being  fitted.  The  ratio  ^ /<3^ 
is  called  the  weighting  factor.  It  can  be  shown  that,  if  the  errors  have 
a  Gaussian  distribution,  by  choosing  the  (^i\    in  this  way  we  get  the 
most  probable  value  of  the  required  D.  We  thus  have  to  solve  the  set  of 
six  linear  algebraic  equations: 

for  the  coefficients  tZ^  :    and  in  particular  the  coefficient  ^60   • 
These  equations  can  be  written  in  the  form 

^  •  -i   -=   Z 
where  !^     is  a  6  x  6  matrix  s.nd  £_    and  ^  are  sis  component  vectors  _, 
the  £-     consisting  of  the  sis  unlsnowns  (Xi.    and  ^    and  _?  being  deter- 
mined by  the  observations. 

By  triangularizing  the  matrix  iJ  we  can  reduce  the  equations  to 
the  form 

where  U       is  a  lower  triangular  matris.   If,  now  ^^^      ,    ^00     and  ?^ 
are  the  leading  elements  of  the  three  matrices  _y   ,  ^  and  ^  ,  respec- 


tively, we  have  directly  that 
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To  determine  the  weighting  factor  ^/^w   i*  is  necessary  to  know 
the  values  of  the  two  standard  deviations.  Owing  to  the  numerous  ways 
errors  can  enter  into  the  observations  and  also  because  of  the  use  of  the 
geostrophic  assumption  it  was  considered  best  to  obtain  its  value  by  the 
experiment  described  below. 

A  series  of  objective  analyses  was  computed  for  the  500  mb  surface 
for  November  25,  1950,  at  0300  GCTo  The  grid  used  is  shown  in  Figure  1. 
Each  ar^^lysis  was  made  with  a  different  weighting  factor.  The  results 
of  each  analysis  were  then  compared  with  the  mean  of  three  subjective 
analyses  for  the  same  map,  prepared  by  three  different  anslysts.  The  re- 
sults are  in  the  form  of  D  values  for  222  points  and  Y^_J)   ,  the  finite 
difference  Laplacian  of  these  D  values,  for  152  points.  The  points  used 
were  selected  prior  to  the  experiment  as  being  those  where  sufficient 
data  were  available  to  permit  a  fair  comparison.  The  means  of  the  stand- 
ard deviations  of  D  and   IP  __[)     of  the  objective  analyses  vs.  the  sub- 
jective analyses  are  shown  on  page  III-22.  It  is  evident  from  that  figure 
that  the  objective  analysis  has  the  best  agi'eement  with  the  subjective 
analyses,  so  far  as  height  is  concerned,  when  wind  receives  a  relatively 
small  weight.  However,  a  heavier  weighting  on  wind  produces  a  better 
agreement  when  W  J^     is  considered.   It  was  believed  best  to  use  a 
weighting  factor  giving  good  results  for  vorticity  calculations,  in  view 
of  the  fact  that  the  objective  analyses  are  principally  for  use  in  nu- 
merical prediction.   It  may  be  noted  that  in  numerical  prediction  one  is 
usually  more  interested  in  expressions  such  as  v    ~1^      and  other 
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111-21.  Location  of  19  x  19  grid  of  points  used  in 
computation. 
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III-22.  Standard  deviations  between  objective  and 
subjective  analyses  for  various  weighting 
factors. 
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derivatives  of  D  rather  than  the  absolute  value  of  D.  Consequently,  a 
weighting  factor  corresponding  to  a  ratio  of  50  ft.  height  error  to  ten 
knots  wind  error  was  selected  for  routine  use  in  upper  air  analyses. 

The  above  t3rpe  of  test  was  not  performed  for  other  levels  but  a 
smaller  weighting  factor,  corresponding  to  30  ft.  height  error  for  a  ten 
know  wind  error,  was  selected  for  use  in  1000  mb  analysis.  This  decision 
was  made  because  of  the  existence  of  greater  errors  in  the  geostrophic 
assumption  near  the  ground.  At  high  levels  wind  observations  are  more 
accurate  with  regard  to  direction  than  to  speed.  This  would  suggest  that 
in  future  analysis  schemes  we  should  use  separate  weighting  factors  for 
the  two  parts  of  a  wind  observation. 

An  analyst,  when  drawing  an  isobar,  looks  at  all  the  wind  and  pres- 
sure data  in  an  area  surrounding  the  small  region  in  question.  In  a  simi- 
lar manner,  the  objective  analysis  program  fits  the  quadratic  surface  as 
closely  as  possible  to  all  wind  and  pressure  data  reported  within  a  pre- 
scribed area  surrounding  the  grid  point.  The  limits  of  the  area  are 
defined  by  the  relation 

|xl   +   |y|   =   lal   , 
where  x  and  y  are  the  coordinates  of  each  piece  of  data  with  respect  to 
the  grid  point  at  which  D  is  to  be  computed.  This  defines  a  square 
(standing  on  one  corner)  with  the  grid  at  the  center.  The  search  then 
consists  of  locating  and  assembling  in  the  memory  of  the  computer  all 
the  data  reported  within  this  area.  These  data  are  then  used  in  the 
fitting  of  the  quadratic  surface  and  the  computation  of  the  D  value  at 
the  grid  point.  The  search  is  made  for  each  grid  point  using  a  given 
sized  area.  When  sufficient  data  are  assembled  for  the  computation, 
the  D  value  is  computed.   If  the  data  are  inadequate  in  quantity,  the 
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point  is  passed  and  a  search  ig  made  around  the  next  point.  In  this  way 
the  computation  proceeds  until  all  points  in  the  grid  have  been  considered. 
This  leaves  some  points  computed  and  some  not  computed. 

The  next  step  is  to  go  back  to  the  points  which  were  skipped  over 
because  of  lack  of  data,  and  to  try  to  assemble  more  data  for  the  computa- 
tion. This  is  done  by  assembling  all  reported  data  plus  all  previously 
computed  data  which  lie  in  the  search  area.  The  search  area  may  be  en- 
larged also.  Thus,  some  of  the  points  which  were  skipped  the  first  time 
the  grid  was  scanned  can  now  be  computed.   (Points  previously  computed  are 
now  skipped  over  and  left  unchanged.)  Some  uncomputed  points  are  still 
left  uncomputed,  due  to  lack  of  data.  After  the  machine  has  proceeded 
through  the  grid  a  second  time,  the  search  area  is  enlarged  again  and  the 
above  process  is  repeated.  The  search  routine  which  is  currently  in  use 
has  evolved  through  many  trials  and  is  summarized  in  Table  I. 


Passage  through 

Grid 

Size  of  Search 

area 

Data  included  in 

(Length  of  side  of 

search 

of 

square,  =  a 

J2) 

1st 

1000  km 

Reported  data  only 

2nd 

1000  km 

^ 

3rd 

1200  km 

l^th 

lUOO  km 

Reported  data  plus 

5th 

l800  km 

/  previously  computed 

6th 

l800  km 

data. 

7th 

l800  km 

J 

Table  I 
The  sudden  increase  in  the  size  of  the  search  area  from  the  l^th  to 
the  5th  search  is  due  to  the  fact  that  by  the  end  of  the  ll-th  search  the 
only  areas  remaining  unanalyzed  are  those  with  practically  no  data  at  all. 
Searches  5  through  7  are  designed  to  fill  in  these  areas  with  values  which 
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are  not  unreasonable. 

Strictly  speaking,  one  needs  six  pieces  of  meteorological  informa- 
tion (each  wind  gives  two  pieces  of  infojnnation  and  each  height  value 
gives  one)  to  fit  a  quadratic  surface  in  the  manner  described  above. 
However,  in  view  of  the  fact  that  the  data  contain  errors,  a  surplus  of 
data  is  highly  desirable  in  order  that  some  smoothing  can  be  done.   An- 
other reason  for  requiring  suprlus  information  is  that,  with  some  combina- 
tions of  data,  a  nearly  singular  matrix  y  is  encountered  in  the  computa- 
tion of  the  D  value.  A  truly  singular  matrix  has  not  jet   been  encountered 
in  any  of  the  analysis  runs,  but  nearly-singular  matrices  have  been  seen. 
These  give  solutions  for  D  which  are  computed  in  the  non-significant  digits 
of  the  data  and  are  therefore  meaningless.  These  can  never  be  entirely 
avoided,  but  the  probability  of  their  occurrence  can  be  reduced  to  a  neg- 
ligible figure  by  requiring  enough  data.  Experience  to  date  has  indicated 
that  this  desired  result  can  be  obtained  if  a  minimum  of  ten  pieces  of  in- 
formation is  required  as  a  prerequisite  to  the  computation  of  the  D  value. 

The  machine  program  can  be  outlined  in  the  following  manner: 

1.  Part  A  code  is  read  into  the  computer  at  the  start  of  the  analysis. 

2.  Part  A  reads  in  data  cards,  one  by  one,  converts  the  data  from  decimal 
to  binary,  and  stores  it.  The  data  cards  contin  the  index  number  of 
the  station  (latitude  and  longitude  for  ships),  the  wind  direction, 
and  speed,  and  the  height  of  the  pressure  surface  to  be  analyzed 

for  each  reporting  station.  Part  A  then  reads  in  part  B  code. 

3.  Part  B  code  reads  in  the  station  locator  cards,  which  have  the  lati- 
tude and  longitude  for  the  corresponding  index  number  of  each  station 
which  can  possibly  report.   It  matches  this  information  with  the  data 
already  in  the  memory  and  converts  the  index  numibera  of  the  data  to 
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latitude  and  longitude.  All  data  are  next  converted  to  the  coordinate 
system  of  the  grid.  Part  B  code  then  reads  the  code  for  parts  C,  D, 
and  E  into  the  machine . 
k.     Part  C  code  is  the  search  routine.  Whenever  enough  data,  are  found 
to  enable  the  computation  of  a  point,  part  C  code  transfers  the  con- 
trol to  part  D. 

5.  Part  D  is  the  matrix  code,  i.e.  the  section  vhich  triangularizes  the 
matrix  y,  etc.  It  computes  the  D  value  at  the  grid  point,  using  the 
data  obtained  in  the  search.   It  then  transfers  the  control  back  to  the 
search  code. 

6.  Part  E  is  the  code  which  converts  the  results  of  the  analysis  from 
binary  to  decimal  and  punches  the  results  on  cards,  which  can  then 
be  tabulated. 

It  can  be  seen  that  the  program  is  entirely  automatic.  The  process 
of  making  the  objective  anal3rsi3  thus  proceeds  in  the  following  way: 

1)  The  reported  meteorological  data  are  punched  on  cards. 

2)  The  cards  are  run  through  a  printer  to  produce  a  listing  of 
the  data.  This  listing  is  then  inspected  by  an  analyst  with  the  aid  of 
a  plotted  map  of  the  data  for  the  detection  of  any  serious  errors  in  the 
data.  This  step  is  unfortunately  necessary  because  of  communication  er- 
rors. However,  it  need  not  delay  the  procedure  unduly,  since  the  maps  can 
be  plotted  at  the  same  time  as  the  cards  are  being  punched.   It  is  not 
difficult  to  postulate  an  addition  to  the  machine  program  which  would  check 
the  data  and  reject  those  containing  serious  errors. 

3)  After  any  serious  errors  detected  in  the  data  are  corrected, 
the  data  cards  are  added  to  the  code  cards,  the  resulting  deck  is  placed 
in  the  computer,  and  the  computation  is  started. 


III-27. 

k)     The  machine  computation  then  proceeds  entirely  automaticallyo 
At  the  present,  slightly  under  1  l/2  hours  are  required  for  the  computa- 
tion of  the  "D"  values  on  one  pressure  surface  for  a  19  x  19  grid  of 
points.   It  is  known  how  to  rewrite  the  code  with  one  slight  internal 
change  which  would  reduce  the  above  computation  time  to  under  one  hour. 
Eventual  reduction  to  a  time  of  less  than  35  minutes  is  visualized  for 
the  IAS  computer. 

5)  The  output  cards  from  part  E  are  usually  printed  to  enable  an 
inspection  of  the  analysis.  These  cards  are  also  used  directly  as  input 
information  for  the  forecast  problem.  The  printed  results  at  each  grid 
point  consist  of: 

a)  the  D  value  obtained  from  the  anslysis,  in  tens  of  feet, 

b)  a  two  digit  number  which  tells  the  number  of  pieces  of 
meteorological  information  used  in  the  computation  of  the  D  value,  and 

c)  a  digit  indicating  at  which  passage  through  the  grid 
the  computation  was  made. 

The  subjective  analyses  available  for  comparison  with  the  objective 
analysis  of  500  mb,  25  November  1950,  0300  GOT,  were  prepared  by  three 
different  individuals,  and  will  be  referred  to  as  A,  B,  and  C  Analyses 
A  and  B  were  prepared  very  carefully  by  experienced  analysts,  who  took 
all  the  time  they  needed  for  their  completion.  Analysis  C  was  made  hur- 
riedly under  operational  conditions.  These  analyses  were  then  compared 
with  an  objective  analysis  (O)  for  the  same  map,  using  the  weighting  fac- 
tor corresponding  to  fifty  feet  height  error  and  ten  knots  wind  error. 
The  standard  deviations  of  D  and  W  2>     in  feet  (for  the  same  points 
referred  to  in  the  description  of  the  determination  of  the  weighting 
factor)  of  the  difference  between  the  various  analyses  are  shown  in 
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Table  II. 


Analysis  pairs 

A-B 

A-C 

B-C 

O-A 

O-B 

0-C 

Std.  Dev.  D. 

50 

67 

60 

69 

72 

68 

Std.  Dev.  V^J) 

13 

17 

16 

18 

20 

18 

Table  II 

It  can  be  seen  from  this  table  that  the  objective  analysis  did  not 
agree  with  the  subjective  analyses  quite  as  well  as  the  subjective  analyses 
agreed  with  each  other.  However,  the  relatively  close  agreement  between 
A  and  B  is  something  which  cannot  ordinarily  be  expected  between  subjective 
analyses.   It  should  not  be  supposed  from  these  results  that  the  objective 
analysis  was  inferior  to  the  subjective  analyses,  since  in  this  case  the 
forecast  is  the  ultimate  test  of  the  worth  of  the  analysis. 

Previously  a  2i)--ho\ar  numerical  forecast  had  been  made  on  the  IAS 
computer  for  the  intense  cyclogenesis  of  November  2i<--25,  1950  starting 
from  1500  Z,   Kovember  2k.     The  input  data  consisted  of  the  1000  and  5OO 
nib  charts  prepared  by  analjrst  B  (above) .  The  same  forecast  was  then  made 
from  objective  analjraes.  The  two  resulting  forecasts,  expressed  in  the 
form  of  2U-hour  forecast  changes  of  height  of  the  1000  mb  surface,  were 
then  correlated  with  each  other  yielding  a  correlation  coefficient  of 
0.96.  Both  forecasts  computed  the  cyclogenesis  which  developed  with  fair 
accuracy,  neither  forecast  being  observably  superior  to  the  other.  From 
this  experiment  it  was  concluded  that  where  a  good  data  coverage  is  found, 
e.g.  over  Korth  America,  two-dimensional  objective  analyses,  prepared  for 
two  different  levels,  are  satisfactory  for  use  in  numerical  prediction. 

When  three  pressure  surfaces  in  the  vertical  are  to  be  analysed 
for  use  in  a  forecast,  a  new  problem  is  encountered.  This  is  the  problem 
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of  representing  accurately  the  horizontal  variation  of  vertical  stability. 
This  variation  is  the  additional  feature  treated  in  the  three -parameter 
model,  as  compared  with  the  two-parameter  model.   In  a  test  of  the  repre- 
sentation of  stability  by  the  above-described  objective  analysis  scheme, 
maps  were  prepared  for  1500  Z,  November  2k,   I95O,  of  a  quantity  directly 
related  to  vertical  stability,  i.e. 


A^T)  =3)„„.-^Xo„^  ">.<>< 


where  the  subscripts  indicate  the  pressure  sirrface  at  which  the  D  value 
was  read  from  the  maps.  This  quantity  was  computed  from  data  at  the  three 
pressure  surfaces  at  each  grid  point  from  the  objective  and  also  from  the 
subjective  analyses.  From  the  grid  point  values,  values  of  Z^^  10     were 
interpolated  at  the  locations  of  all  reporting  radiosonde  stations. 
Values  of   A  3^  were  also  calculated  from  the  actual  radiosonde  re- 
ports. These  were  used  as  verification  of  the  values  obtained  from  the 
objective  and  from  the  subjective  analyses.  The  values  of   A.  JD  ob- 
tained by  objective  analysis  when  correlated  with  the  values  from  the 
actual  observations  yielded  a  correlation  coefficient  of  O.9O.  The  cor- 
relation coefficient  between  the  values  obtained  from  subjective  analjrsis 
and  the  reports  was  0.9^«  From  this  it  was  concluded  that  the  two- 
dimensional  scheme  of  objective  analysis  gives  an  adequate  representation 
of  the  horizontal  variation  of  vertical  stability  where  the  data  cover- 
age is  good  and  when  the  pressure  surfaces  analyzed  are  at  least  300  mb 
apart . 

The  next  experiment  consisted  of  using  the  maps  prepared  for 
1500  Z,  2^4-  November  1950,  in  two  forecasts  with  a  three -parameter  model 
--  one  to  be  made  from  the  subjective  analyses  and  the  other  from 
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III-30.  l<-00  mb  initial,  forecast  and  verification 
maps.  The  dashed  line  indicates  the  ob- 
served 2ll-hour  movement  of  the  low  center. 
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III-31.   1000  nib  inioial,  forecast  and  verification 

maps.  The  dashed  line  indicates  the  observed 
pk-hnnr  mnvement  of  the  low  center. 
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the  objective  analyses.  The  results,  at  lj-00  and  1000  mb,  of  these  fore- 
casts are  presented  in  the  figures  on  pages  III-30  and  III-31.  The  700 
mh  results  are  very  similar.  These  figures  show  the  surprising  result 
that  the  forecast  made  from  the  objective  analyses  verified  better  than 
the  forecast  made  from  the  subjective  analyses.  Further  confirmation  of 
this  is  given  by  the  correlation  coefficients  of  forecast  change  vs.  ob- 
served change  at  400  rob,   where  the  difference  between  the  two  forecasts 
was  greatest.  These  are  given  in  Table  III. 


Correlation 
Coefficient 

Forecast  change  vs.  observed 
change — subjective  analysis 

0.63 

Forecast  change  vs.  observed 
change--obJective  analysis 

0.79 

Forecast  change — objective 
analysis  vs.  forecast  change 
— subjective  analysis 

0.82 

Table  III 
The  results  described  above  can  be  attributed  to  the  fact  that 
when  the  objective  analyses  are  made,  the  geostrophic  approximation  is 
forced  on  the  analyses,  thus  giving  better  results  when  a  quasi-geostrophic 
model  is  used  to  make  the  forecast.  This  forcing  is  most  marked  in  the 
vicinity  of  pronounced  troughs  or  lows  at  middle  and  high  levels  in  the 
troposphere,  where  the  speed  and  curvature  of  the  flow  are  strong.   It 
takes  the  form  of  diminishing  the  analyzed  geostrophic  wind  from  the  ac- 
tual geostrophic  wind  in  the  areas  where  the  observed  winds  are  markedly 
sub-geostrophic.  This  is  due  to  the  fact  that  observed  winds  as  well  as 
all  the  heights  of  the  pressure  surface  are  used  in  the  objective  analysis. 
The  smoothing  process  incorporated  into  the  objective  analysis  also 
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contri'buteB  to  the  same  effect,  particularly  when  the  half -wave  length 
of  the  wind  pattern  is  not  much  larger  than  the  dimension  of  the  search 
area.  The  quasi-geostrophic  prediction  model,  using  the  geoatrophic  ap- 
proximation for  measurement  and  advection  of  vorticity,  then  has  a  more 
useful  representation  of  these  quantities  than  it  otherwise  would. 

The  results  obtained  above  lead  to  the  conclusion  that  a  simple 
two-dimensional  scheme  of  objective  analysis  is  feasible  and  satisfactory 
for  the  purpose  of  transforming  raw  data  into  input  data  for  a  two-  or 
three-parameter  numerical  prediction  model.  In  fact,  such  a  scheme  ap- 
pears to  be  better  suited  for  use  in  a  quasi-geostrophic  prediction  model 
than  information  obtained  from  subjective  analysis.  The  principal  dif- 
ficulty encountered  so  far  is  introduced  by  incomplete  and  erroneous 
meteorological  observations.  The  detection  and  correction  of  the  er- 
rors at  present  must  be  done  subjectively,  i.e.,  a  skilled  analyst  must 
inspect  the  data  very  carefully  for  omissions  and  errors.  One  could  argue 
that  this  process  practically  amounts  to  a  subjective  analysis.  It  should 
be  pointed  out  that  in  the  preparation  of  input  data  for  a  forecast  from 
a  subjective  analysis,  the  analysis  itself  is  only  the  beginning  of  the 
process.  The  reading,  recording,  and  checking  of  interpolated  informa- 
tion at  the  grid  points  is  the  most  tedious  and  lengthy  part  of  the 
process.  All  this  is  unnecessary  when  an  objective  analysis  method  is 
used.  It  is  clear,  however,  that  the  reduction  and  eventual  elimina- 
tion of  errors  of  computation  and  transmission  of  meteorological  data  is 
a  subject  which  must  receive  high  priority  in  future  research  and  develop- 
ment before  a  truly  objective  analysis  is  possible. 
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H.  Wind  driven  vater  waves. 

Several  investigations  have  been  undertaken  in  the  past  to  deter- 
mine the  character  of  the  inertial  oscillations  initiated  by  an  unbal- 
anced current  in  the  oceans.  These  investigations  have  been  restricted 
to  either  the  final  equilibrium  state  or  to  a  system  in  which  the  momentum 
causing  the  unbalanced  current  has  been  added  impulsively.  The  present 
investigation  is  concerned  with  the  effect  on  the  resultant  motion  of  the 
length  of  time  during  which  the  momentum  is  added. 

A  simple  model  is  constructed  as  follows:  A  rotating  basin  of  homo- 
geneous fluid  of  infinite  horizontal  extent  and  of  finite  depth  is  initi- 
ally at  rest  relative  to  the  rotating  coordinate  system.  At  time  t'  =  0, 
a  wind-stress  is  applied  to  an  infinite  strip  of  the  surface  and  continues 

to  act  until  time  t'  =  t  '.  \fha.t   are  the  velocity,  mass,  and  energy  dis- 

o 

tributions  at  time  t? 

If  the  x'  and  y'  axes  are  directed  along  the  horizontal  and  the  z' 
axis  is  taken  positive  upward,  the  above  problem  reduces  to  the  following 
system  of  equations i 

^jf     ^   /^    =   -  ^5  3>  ^ 

where  u  and  v  are  the  (vertically)  integrated  mass  transports  in  the  x' 
and  y'  directions,  respectively, 

n  is  the  height  of  the  water  above  the  equilibrium  surface, 

D  is  the  initial  depth, 

f  denotes  the  (constant)  Coriolis  parameter, 
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•JT  is  the  wind-atreas  applied  at  the  surface, 

g  is  the  gravitation  acceleration,  and 
t'  ia  the  time. 

The  above  set  of  equations  can  be  solved  analytically  when  the 
wind-stress,  "C      ,  is  a  sufficiently  simple  function.  However,  the 
task  of  computing  the  results  by  hand  is  formidable  and  the  electronic 
computer  was  used.  The  computational  problem  is  considerably  simplified 
by  integrating  the  original  three  differential  equations  directly.  In 
order  to  perform  the  above  integration  the  equations  were  first  made  non- 
dimensional. 

The  wind-stress,  "Z^  ,  is  of  the  form 
r:,.  =  W        t'  <  t^',   y  <  a  , 
r^  =  0        t'  >  t^',   all  y, 
r^  =  0        t'  <  t^',   y  >  a. 
The  non-dimensional  equations  are  then  given  by  defining 

li^  _  \/  :^  r- 

l^   ^  0  -  -^ 

T   =         :^     ;     h)   <     '  ^      0<    ±   <     to 

r-  =     o      ^    I::))  >  I     ^    o<  -^  <  -to 

T  -       O  .      i=  >  ir^ 
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In  centered  finite  difference  form  the  equations  are 


A+H,    =   -Ai    A^v"" 


where  the  superscript,  0  ,   refers  to  the  time  index  and  the  subscript, 
±,   refers  to  the  space  point.  Thus 

The  equations  can  now  be  written 

An  inspection  of  the  equations  reveals  the  fact  that  ^^ti  does 
not  appear  in  the  equation  for  A^^  .  Thus,  the  space  points  can  be 
staggered  as  shown  below 

— ?^ — >< — ^ — f — n >  y 

In  this  case,  if  the  U  and  V  values  are  stored  together  with  the 
following  H  value,  the  computation  time  can  be  reduced  by  a  factor  of  2. 
The  equations  then  are  changed  by  the  definition 
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The  numerical  values  of  the  parameters  are 

2  -h       -1 

a  =  100  km,  g  =  lOm/sec  ,  D  =  i^■000  m,  f  =  10  sec 

^5  =  .5  (=  50  km  in  dimensional  form),   ^      =  l/^OO, 

AT  =   .0238  (=  1/22  pendulum  days  in  dimensional  form). 

These  values  are  chosen  so  as  to  satisfy  computational  stability 

criteria.  Because  of  the  wide  range  of  values  of  the  variables,  the 

5  —  -5  — 

equations  were  scaled  by  defining  V  =  2  -^  V  for  all  i,  and  U  =  2   U  for 

i  >  0  (only  one  point  i  =  0  is  taken  in  the  region  within  which  the  wind- 
stress  is  acting) .  The  two  sets  of  equations  in  scaled  form  are  thus 


U^  =     U^  +  J  AT-  1/ 


^^-  _-  VZ-'-^ATUf-^^fT^  .. 


As  ^^ 
A' AS 
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Because  of  the  rapidity  with  which  the  waves  move  to  the  outlying  regions, 
some  means  of  confirning  the  computation  to  the  region  of  interest 
(O  <  i  <  100)  is  necessary.  This  was  accomplished  by  considering  two 
syBtems  of  the  above  tjrpe  with  opposite  boundary  conditions  at  the  ex- 
tremity. The  symmetry  of  the  system  of  interest  about  the  y  axis  enables 
one  to  specify  boundary  conditions  which  cancel  when  added.  Because  of 
the  linear  character  of  the  equations  the  solutions  could  be  superposed. 

The  two  problems  treated  are  shown  schematically  below  (looking 
down  on  the  surface  of  the  ocean) 

ITI  I  It 

When  these  two  problems  are  added  (note  that  the  system  on  the  right  in 
each  case  can  be  accounted  for  by  the  specification  of  a  boundary  condition 
along  the  dotted  line  half-way  between  the  two  systems)  the  effect  of  the 
right  half  cancels  out  and  one  is  left  with  twice  the  values  of  the  left- 
hand  system. 

Computations  for  2^  pendulum  hours  were  carried  through  for  three 
cases.  The  cases  are  distinguished  by  letting  t  take  on  different  values. 
Case  I     t  =  0   corresponding  to  0  pendulum  hours 
Case  II    t  =  6   corresponding  to  6  pendulum  hours 
Case  III   t  =  12  corresponding  to  12  pendulum  hours. 
The  results  indicate  a  pronounced  dependence  of  the  amplitude  and 
energy  of  the  oscillations  on  the  time  of  duration  of  the  wind-stress. 
It  is  shown  that  if  the  wind  acts  for  12  hours,  the  oscillations  are  con- 
siderably smaller  in  amplitude  than  those  in  which  the  wind  is  added  im- 
pulsively. 
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The  above  problem  vas  also  considered  for  an  ocean  consisting  of 
two  homogeneous  layers  of  different  density  with  the  lower  layer  contin- 
ually adjusting  in  such  a  way  as  to  offset  any  pressure  gradient  imposed 
by  the  surface  of  the  ocean  (such  an  assumption  is  justified  when  the 
lower  layer  is  much  deeper  than  the  upper) .  The  results  in  this  case 
show  that  the  amount  of  energy  absorbed  by  the  oscillations  is  very  much 
dependent  on  the  length  of  time  of  wind-stress  duration.  In  fact,  a  wind 
acting  for  a  period  of  12  hoiors  sets  up  oscillations  which  have  an  energy 
which  is  about  3^  of  that  absorbed  by  oscillations  caused  by  an  impulsive 
addition  of  momentum. 
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I.  Use  of  the  primitive  equations  in  foreeaetlng. 

Because  of  the  limitations  Imposed  on  the  accuracy  of  the  fore- 
cast by  the  use  of  the  geostrophic  approximation,  several  attempts  have 
been  made  recently  to  forecast  the  pressure  and  velocity  fields  by  means 
of  the  primitive  equations .  By  thus  allowing  gravity  waves  to  enter  into 
the  solution,  one  encounters  severe  difficulties.  Specifically,  the  ef- 
fects of  the  boundary  conditions  are  felt  throughout  the  domain  within 
three  hours.  As  the  boundary  conditions  are  necessarily  of  a  fictitious 
nature,  the  large-scale  meteorological  flow  is  obscured  and  the  forecasts 
are  practically  worthless. 

It  has  been  suggested  that  a  simple  model  be  prescribed  so  that 
stability  of  boundary  conditions  may  be  easily  studied.  The  following 
is  a  report  of  one  such  attempt. 

The  domain  Is  defined  by  a  rotating  cylinder  extending  from  y  =  0 
to  y  =  T  along  the  generator  and  from  x  =  0  to  x  =  X  circumferentially. 
The  points  x  =  0  and  x  =  X  coincide,  i.e.,  the  domain  is  finite  and  the 
motions  are  periodic  circumferentially.  The  boundaries  at  y  =  0  and  y  = 
Y  are  solid  boundaries. 

The  primitive  equations  for  the  barotropic  model  are: 


Tt   "   ^^^  -^  V  V 


III-ln. 


■where  u  and  v  are  the   velocities  along  the  x  and  y  directions,  respectively, 
m      is  the  geopotential  height  and  f,  the  Coriolis  parameter,  is  constant. 
In  finite  difference  form,  the  scaled  equations  are 

vhere  the  reference  point  has  subscripts  i,  j  (O  <  i  <  I,  0  <  j  <  J) 

y/^  ,  /xT  >   A  5  '  "B  '  ^  '  ®'^®  scaling  constants.  Superscripts  refer 
to  time  index.  ^C/^-       ,  '^ij      >     rcj       s^e  scaled  quantities  (modulus  <  1)  • 

For  interior  points  of  the  domain  equation  s  {k) ,  (5) ,  (6)  suffice 
for  evaluating  the  velocities  and  the  pressure  for  future  times.  The  same 
holds  true  for  boundary  points  at  x  =  0,  X  since  these  point  may  also  be 
considered  "interior"  because  of  the  periodicity. 

At  the  boundaries  y  =  0,  y  =  Y,  the  normal  derivatives  cannot  be 
given  by  centered  space  differences.  It  can  be  seen  from  equation  {k) , 
however,  that  no  difficulty  arises  in  evaluating  ^,j  since  each 

derivative  normal  to  the  boundary  is  multiplied  by  ^'   which  vanishes 
at  these  points.  The  term  y(/^        is,  of  course,  zero  at  all  times  at 
the  solid  boundaries.  One  is  left,  finally,  with  the  problem  of  evaluat- 
ing  ©.   at  y  =  0,  Y.  Two  methods  have  been  considered. 
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(7) 


Case  I.  Equation  (5)  written  for  C  +  1   reduces  to 

where  +  depends  on  i  =  I  or  0,  respectively.  Thus  <p  ^j      can  be  evalu- 
ated easily  and,  in  fact,  this  method  indicates  definite  stability  for 
a  forecast  of  2l<-  hours. 

Case  II.  In  eqaution  (6)  the  only  unknown  term  is  A<j  i^  •  If  > 

however,  we  choose  as  Z^uiT  ,  twice  the  one-sided  derivative  of  ^</  ^j 

0  r+i 

at  the  boundary,  then  <p^  ■         can  be  evaluated  directly.  Case  II  results 

are  not  as  good  as  those  of  Case  I. 

A  third  possibility  involves  taking  one-sided  space  and  time  deriva- 
tives for  all  quantities  in  equation  (6) .  This  method  has  not  as  yet  been 
carried  through. 
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